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0\ : ABSTRACT 

We present analytical approximations for calculating the scattering, absorption and 
escape of non-ionizing photons from a spherically symmetric two-phase clumpy medium, 
Q\ ! with either a central point source of isotropic radiation, a uniform distribution of isotropic 

emitters, or uniformly illuminated by external sources. The analytical approximations 
are based on the mega-grains model of two-phase clumpy media, as proposed by Hobson 
& Padman, combined with escape and absorption probability formulae for homogeneous 
media. The accuracy of the approximations is examined by comparison with three- 
dimensional Monte Carlo simulations of radiative transfer, including multiple scattering. 
Our studies show that the combined mega-grains and escape/absorption probability 
formulae provide a good approximation of the escaping and absorbed radiation fractions 
for a wide range of parameters characterizing the dumpiness and optical properties of 
5-h 1 the medium. 

A realistic test of the analytic approximations is performed by modeling the absorp- 
tion of a starlike source of radiation by interstellar dust in a clumpy medium, and by 
calculating the resulting equilibrium dust temperatures and infrared emission spectrum 
of both the clumps and the interclump medium. In particular, we find that the tem- 
perature of dust in clumps is lower than in the interclump medium if the clumps are 
optically thick at wavelengths where most of the absorption occurs. Comparison with 
Monte Carlo simulations of radiative transfer in the same environment shows that the 
analytic model yields a good approximation of dust temperatures and the emerging UV 
to FIR spectrum of radiation for all three types of source distributions mentioned above. 

Our analytical model provides a numerically expedient way to estimate radiative 
transfer in a variety of interstellar conditions and can be applied to a wide range of 
astrophysical environments, from clumpy star forming regions to starburst galaxies. 



Of 



Raytheon ITSS, varosi@gsfc.nasa.gov 



Varosi & Dwek: Escape and Absorption of Radiation in Clumpy Dusty Environments 



Contents 



7 MODELING THE ABSORPTION OF 



STELLAR RADIATION BY DUST AND 



INTRODUCTION 



THE INFRARED EMISSION 



MONTE CARLO SIMULATION Oi 



7.1 Absorbed Luminosities and the 



37 



38 



RADIATIVE TRANSFER 



2.1 General Simulation Methods 



5 



Distribution ol Dust Temperatures 



7.2 The Emerging UV-FIR SEP 43 



2.2 Radiative Transfer in Two-Phase Clumpy 



7.3 Exploring Parameter Space with the 



Media 



MGEP Model 



44 



3 ESCAPE AND ABSORPTION IN 



HOMOGENEOUS MEDIA 



8 SUMMARY AND CONCLUSIONS] 46 



10 



|3 . 1 Central Isotropic Point Source 10 



A FINAL FLUX APPORTIONMENT IN 



MONTE CARLO 



47 



£L3 — Finite Number of Uniformly Distributed! 


B RANDOM CLUMPS IN A TWO PHASE 


Point Sources 


J 13 MEDIUM 




3.4 Uniformly Illuminating External Source 14 


C ESCAPE AND INTERACTION PROB- 



49 



4 ' ESCAPE — A-N-B — ABSORPTION — IN 



CLUMPY MEDIA 



ABILITIES FOR A SPHERE 



49 



4.1 The Mega- Grains Approximation 



4.1.1 Effective Optical Depth 



4.1.2 Effective Scattering Albedo 



16 

16 
17 
17 



|C.l No Scattering! 49 

|C.2 Including Scattering 51 



D ANALYSIS OF CLUMP OVERLAPS 52 



4.1.3 Effective Scattering Asymmc- 



try Parameter] 

^^^^n^^xtendc^^ega-Grains Approxi- 



19" 



E DISTRIBUTION OF DUST TEMPER- 



ATURES AROUND A POINT SOURCE 



mation 



OF RADIATION 



53 



4.3 Comparison with Earlier Theory 
[4.4 Escape and Absorption Probabilities 



20 
2G 



for Clumpy Media 27 



[4.5 The Fractions Absorbed in Each Phase 



of a Clumpy Medium 28 



5 SUMMARY OF EQUATIONS 



5.1 Escape and Absorption Probabilities 



5.2 Mega-Grains Approximation 



5.3 Fractions Absorbed in Each Phase 



5.4 The Case of a < 1 



29 

29 
29 
30 
30 



S COMPARISON OF ANALYTIC APPROX 



IMATIONS WITH MONTE CARLO SIM 



ULATIONS 



30 

3.1 Dependence on Optical Depth 31 



3.2 Dependence on Filling Factor 33 

3.3 Dependence on Density Ratio .... 35 



3.4 Body Centered Cubic Lattice of Clumps 35 



Varosi & Dwek: Escape and Absorption of Radiation in Clumpy Dusty Environments 



3 



1. INTRODUCTION 

Radiative transfer plays an important role in the 
spectral appearance of almost all astrophysical sys- 
tems ranging from isolated star forming regions to 
protogalactic and galactic systems. At wavelengths 
longer than the Lyman limit, absorption and scatter- 
ing by dust in the intervening medium is the major 
fac |6r affecting the transfer of radiation. Most of the 
transfer of radiation occurs in the interstellar medium 
(ISM) of the host system of the emitting sources. For 
simplicity, models of radiative transfer often assume 
a homogeneous distribution of dust and gas, although 
in reality the structure of the ISM is observed to be 
significantly more complex. 



' n our Galaxy, the ISM is known to be composed 



of at least five phases that are in approximate pres- 
sure equilibrium: cold dense molecular clouds, cold 
diffuse clouds, warm diffuse clouds, H II regions, and 
hot low density cavities created by supernova rem- 



nants ( jBpitzcr 1978| , |Cox 1994 |McKee 1994 |Knapp 
1995). The ISM is observed to have clumpy struc- 
ture even at very small scales ( Marscher et al. 1993| , 
stutzki & Gusten 1990). Data from CO line emis- 



sion indicate a distribution of molecular clouds having 



a power-law cloud mass spectrum ( 


Sanders, Scovillc, 


and Solomon 1985 


)■ 


Dickey & Garwood (1989) find 



the same cloud mass spectrum based on 21 cm emis- 
sioih line data. Most likely the ISM has a spectrum 



of densities and temperatures with correlated multi- 
scale spatial structure, as evidenced by sky surveys 
such as the Infrared Astronomical Satellite (IRAS) 
and 21 cm surveys. Analysis of IRAS 100/^m sky 
flux indicates that the diffuse H I clouds have a frac- 



tal distribution (Bazcll fc Desert 1988, Waller et al 



1997 ). Analysis of CO column densities and emission 
line profiles further suggests that a fractal distribution 
of matter applies to the molecular component of the 



ISM as well (Falgarone 1995, Elmegreen & Falgaronc 



1996). The complex, possibly fractal, distribution of 



gas and dust is supported by theoretical arguments 



as well ( 


Pfenniger & Combes 1994 




Rosen & Breg- 


nan 1994 Norman & Ferrara 1996 




Elmegreen 1997 


)• 



Thus the ISM is clearly inhomogeneous and simula- 
tions of radiative transfer should account for this fact. 

The simplest model of an inhomogeneous medium 
is that consisting of two phases: dense clumps em- 
bedded in a less dense interclump medium (ICM). 



Natta & Panagia (1984) developed simple analytic ap- 



mcdia with no scattering and an empty ICM. Radia- 
tive transfer with isotropic scattering in a two-phase 
clumpy medium (non-empty ICM) with plane-parallel 
geometry was investigated by Boissc (1990)| . He used 
a Markov process model of the medium to develop 
analytical approximations for the intensity of radia- 
tion. Comparison with 3D Monte Carlo simulations 



verified the accuracy of the approximations. Hobson 



& Scheuer (1993) developed analytic solutions of ra- 
diative transfer with isotropic scattering in A— phase 
clumpy media using a Markov process model, which 
they compared with Monte Carlo simulations for the 
cases of two or three phases. Analytic approxima- 
tions for radiative trans fer in two-phase clu mpy media 
were also develo ped by Neufcld (1991)| and Hobson & 



Padman (1993), and we discuss and utilize them in 



later sections of this paper. Witt fc Gordon (1996) 
performed extensive Monte Carlo simulations of the 
transfer of radiation from a central point source in a 
spherical two-phase clumpy medium when the scat- 
tering is not isotropic, more typical of UV photons 
scattered by dust. Their simulations showed that the 
nonlinear variation of effective optical depth and ef- 
fective albedo with respect to parameters characteriz- 
ing the clumpy medium could lead to erroneous esti- 
mates of the dust albedo and opacity if one assumes a 
homogeneous medium when modeling what may ac- 
tually be a clumpy medium. Gordon, Calzetti, and| 



Witt (1997)| applied the Monte Carlo model of Witt 
& Gordon to the study of dust in starburst galax- 
ies, concluding that a shell of clumpy dust around a 
starburst provides the best model of spectral data. 
Radiative transfer in a clumpy environment was also 
investigated by Wolf, Fischer & Pfau (1998), with re- 



sults similar to Witt & Gordon. 

All the simulations and studies mentioned above 
verify the general expectation that a medium is more 
transparent when it is clumpy. This can be demon- 
strated in the following generic example, which then 
leads to the concept of effective optical depth in an 
inhomogeneous medium. Consider N randomly cho- 
sen parallel and equal length lines of sight of through 
an inhomogeneous medium acting as a foreground 
screen. Let Tj be the optical depth of the i-th. line of 
sight, defined as the product of column density and 
cross-section of the dust, so that e~ Ti is the trans- 
mission. If N is large (e.g. the number of photons 
per second emitted by a galaxy) then we can com- 
pute the average transmission of all the lines of sight, 



proximations for the effective optical depth of clumpy 
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wh ch defines the effective optical depth: 



1997). For these types of modeling efforts, easily com- 



T eff 



= - ln (^E ex p(~ T *)) 



(i) 



On the other hand, the average of all the optical 
depths equals the optical depth of the homogeneous 
medium with equal dust mass: 



The 



N ^ 



(2) 



Using standard calculus one can prove the inequality 



(3) 



which states that the average transmission of the 
inhomogeneous medium is greater than that of the 
equivalent homogeneous medium, allowing relatively 
more photons to escape. The expression can be an 
equality only when n = Thom for all lines of sight, 
i.e. only if the medium is homogeneous. Applying 
the negative natural logarithm to the inequality (^) 
and using the definitions (Q) and (||) gives 



r eff < T ho 



(4) 



so the effective optical depth of an inhomogeneous 
medium is less than that of a homogeneous medium 
with equal mass of dust. When the dust also scat- 
ters photons, the above inequality can be consid- 
ered to apply approximately to the scattered photons, 
so we again expect generically a greater transmis- 
sion through inhomogeneous media than the equiv- 
alent homogeneous medium. In this paper we study 
the dependence of the nonlinear relationship between 
r e f f and Thom on parameters characterizing a clumpy 
medium, and also give approximations for how ab- 
sorption of photons is apportioned in each phase of 
the medium. 

One consequence of the fact that r e ff < r^ om 
for an inhomogeneous medium is that a large mass 
iuot could bo concealed in donoo rogiono oven 



of 

tho ugh oboorvationo of low extinction infer a email 



mass of dust. Such dense clumps of dust could have 
a lower temperature than if the dust were distributed 
homogeneously, providing an alternative explanation 
of observations that show an inconsistency between 
dust temperature, absorbed luminosity, and known 
luminosity sources, as found in recent infrared spec- 
tral observations of the Galactic Center ( phan ct al 



putable analytic approximations of radiative transfer 
in clumpy media are desirable since they would en- 
able the rapid exploration of the effects of varying 
the parameters in an astrophysical model. The more 
time consuming Monte Carlo simulations can be used 
to guide the development and test the accuracy of 
analytical approximations. 

We have developed a general Monte Carlo Radia- 
tive Transfer (MCRT) code for simulating radiative 
transfer with multiple scattering in a three dimen- 
sional inhomogeneous medium, utilizing some novel 
techniques. The model is restricted to non-ionizing 
radiation, hence the only absorption and scattering 
considered is due to dust. We repeated many of the 
Monte Carlo simulations of Witt & Gordon and our 
results are in good agreement. The medium in their 
model is composed of cubic clumps randomly located 
on a body centered cubic lattice. We instead con- 
centrate our studies on two-phase media composed of 
spherical clumps because the radiative transfer prop- 
erties of such media are more directly approximated 
by analytic formulas in the mega-grains model of 



Hobson fc Padman (1993), hereafter HP93. Further- 



more, we extend the research of HP93 by comput- 
ing the fraction of photons absorbed in each phase 
of the medium and improving the approximations for 
scattering. The resulting formulae allow for conve- 
nient and rapid exploration of the effects of varying 
the parameters defining the environment on the es- 
cape and absorption of radiation. The cases of a 
central isotropic point source, uniformly distributed 
isotropic emitters, and uniformly illuminating exter- 
nal sources are studied for homogeneous and clumpy 
media, with both MCRT simulations and analytic 
approximations, over a wide range of dust optical 
depths, scattering albedos, and from isotropic to for- 
ward scattering. The validity and accuracy of the 
analytical approximations is tested by detailed com- 
parison with MCRT simulations. 

We have also investigated the transfer o f radiation 
in a fracta l distribution of dust density ( Varosi & 



Dwek 1997) using MCRT, where the fractal construc- 



tion is based on methods given in Elmegreen (1997). 
Briefly, the variation of the escape of radiation in a 
fractal medium as a function of T/ lom is similar to that 
of clumpy media, and we find that the mega-grains 
model can be used to approximate the escaping frac- 
tion of photons if some additional assumptions are 
imposed. However, the distribution of absorbed radi- 
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ation is more complicated than in two-phase clumpy 
media since there is a spectrum of densities in the 
fractal case. Studies and approximations of radiative 
transfer in complex and fractal media will be pre- 



sented elsewhere (Varosi & Dwek 2000). 

The paper is organized as follows. Section 2 de- 
scribes the MCRT code and its application to two- 
phase clumpy media. In we present analytical 
approximations for the fraction of radiation escaping 
from or absorbed by homogeneous media in spherical 
geometry, for each of the three types of sources men- 
tioned above. In §^ we present the mega-grains model 
of two-phase clumpy media, and also some improve- 
ments and extensions. The mega-grains approach 
is combined with the formulae for homogeneous me- 
dia to get approximations for the escaping/ absorbed 
fractions in clumpy media, and also the fractions ab- 
sorbed in clumps and the ICM. Section || gives a sum- 
mary of all the equations in a convenient list. Then 
in §|6|, we demonstrate the validity of the analytic ap- 
proximations by comparison with MCRT simulations. 
In §0we present a realistic test of the analytical model 
by simulating the scattering, escape, and absorption 
of a starlike emission spectrum by a clumpy distribu- 
tion of interstellar dust, and then comparing the pre- 
dicted equilibrium dust temperatures and emerging 
UV-FIR spectrum with the results of MCRT simula- 
tions. The analytical model is also used to explore 
a large region of parameter space characterizing the 
clumpy medium, and we explain why the dust tem- 
perature in clumps is lower than in the ICM. The 
summary and conclusions are presented in §8. 

2. MONTE CARLO SIMULATION OF 
RADIATIVE TRANSFER 

2.1. General Simulation Methods 

The Monte Carlo code we have developed simulates 
coherent scattering and absorption of photons by dust 
with any kind of spatial distribution and composition. 
The spatial distribution of the dust is specified by a 
continuous function p(x,y, z), the mass of dust per 
unit volume. In practice this quantity is defined on a 
discrete three-dimensional grid. For each wavelength 
simulated, the number of photons absorbed by the 
dust in each volume element (voxel) of the 3D grid is 
saved, allowing computation of the dust temperatures 
and resulting infrared emission spectrum. 

The grid resolution is limited only by the available 
computer memory: increasing the number of grid el- 



ements does not significantly affect the computation 
time. This is achieved by employing the Monte Carlo 
method of rejections (Neumann 1951) in selecting the 
random distances each photon travels between inter- 
actions (absorption or scattering) with the dust, as 
described in [Lux fc Koblinger (1995, pp.40-41)|. The 



method proceeds as follows. Let r = (x, y, z) be the 
initial position of the emitted photon, and p max be 
the maximum dust mass density in the direction 1 
the photon is traveling. Assume temporarily that 
the density is uniform and equal to p max along its 
path. In that case the probability that the photon 
will travel a distance s without having an interaction 
is exp(— sK/o max ), where k is the dust mass extinction 
coefficient. Applying the fundamental principle of the 
Monte Carlo method, one can choose a uniformly dis- 
tributed random variable < C < 1 an d set £ to be 
equal to 1 — exp(— sAtp max ), the probability that an in- 
teraction will occur, and then solve for s, the random 
value for the interaction distance, to give: 



ln(l - C) 



(5) 



Note that (/cpmax) -1 is the worst case mean- free-path. 
One then has to play a rejection game in order to de- 
termine if the supposed interaction at the new loca- 
tion, r' = r+sl, is accepted as a real interaction, be- 
cause the density is of course not everywhere equal to 
Pmax- We choose another uniformly distributed ran- 
dom number < i? < 1 and compare it to p(r')/ p max , 
the ratio of the actual density at the event location 
r' to the maximum density along 1. If d < p{r')/p msix 
the supposed interaction is accepted as real. In par- 
ticular, if /o(r')/p max = 1, all values of the random 
variable $ will fall below this ratio, and the interac- 
tion will be accepted as real. Thus p(r')/p max is equal 
to the probability that an interaction is real. If it hap- 
pens that d > p(r')/p max the interaction is rejected 
and called a virtual interaction. After a virtual in- 
teraction the photon is allowed to travel in the same 
direction another random distance selected by eq.(|J), 
and the above steps are repeated until an interaction 
event is accepted as real or the photon escapes. After 
evolving many photons in this manner, the method 
effectively integrates the density along all directions 
of travel while performing the Monte Carlo simula- 
tion of photon interactions. The 3D medium need be 
specified with only a point density function, and the 
method is simpler than performing numerical integra- 
tions along each photon path to determine in which 
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voxel the photon will interact with dust. For most 
simulations one can set p max equal to the maximum 
density of the entire medium, thereby further simpli- 
fying the process. However if /0 max ^> p m i n and the 
high density regions occupy a very small fraction of 
the volume, it is worthwhile to have an algorithm for 
determining the maximum density in a given direc- 
tion, otherwise the method may iterate through many 
virtual interactions before getting to a real interaction 
or escaping, using up more CPU time. 

Our Monte Carlo code follows a large group of ran- 
domly emitted photons simultaneously, typically 10 5 
at a time, through multiple scatterings and to even- 
tual absorption or escape. Each photon is given an 
initial flux weight of unity, and this weight changes as 
the photon interacts with the dust. The simulation 
proceeds via iterations, as follows. First, the group 
of photons travel randomly selected distances (as de- 
scribed above) depending on the dust density encoun- 
tered, and then photons either escape the medium 
(reducing the number in the group) or interact with 
dust. Photons that interact are considered to expe- 
rience both scattering and absorption: they are re- 
duced in flux weight by the dust scattering albedo 
u) = <T SC at I a ext (the ratio of scattering to extinction 
cross-section), and a fraction 1 — lo of the energy of 
each interacting photon is absorbed by the dust at the 
locations of interaction. The path of each scattered 
photon is deflected by a random angle, 9 scat , which is 
distributed according to the Henyey-Greenstcin (HG) 
phase function (Hcnyey & Grccnstcin 1941), fully 
characterized by the parameter g = (cos 9 sca t), the av- 
erage value of the cosine of the deflection angles, also 
called the asymmetry parameter. The random scat- 
tering deflection angles are selected by the method 
given in Witt (1977). The remaining photons hav- 



ing new directions are fed back into the beginning of 
the iteration procedure. This cycle of interaction, ab- 
sorption and scattering is repeated until the net flux 
of photons remaining in the medium is a small frac- 
tion of the flux that has escaped (e.g. less than 5%). 

The flux of photons remaining after the iterations 
are terminated is divided into escaping and absorbed 
fractions according to a new scheme that considers 
the history of the escaping photons during each iter- 
ation and the dust scattering albedo. Let Nk be the 
number of photons remaining after iteration fc, and 
let m be the total number of iterations performed. If 
the remaining fraction after each iteration, given by 
the ratio (3k = Nk/N^-x < 1, is tending toward a 



constant value, then the fraction of the final N rn pho- 
tons that will escape if the iterations were continued 
ad-infmitum is approximately (see Appendix [A]) 



r esc 

J m 



1 — LUj3„ 



(0) 



Thus upon termination after m scatterings we assume 
that a fraction f^ sc of the remaining flux escapes and 
a fraction 1 — /,^ sc is absorbed at the current loca- 
tion. Testing of the formula has shown that even if 
the multiple scattering is terminated after a few iter- 
ations when the flux remaining is 20% of the escaped 
flux, the final escaping fraction of flux obtained after 
applying eq. (|6|) is as accurate as if the iteration was 
continued until the remaining flux is less than 1% of 
the escaped flux, and the total computation time is re- 
duced. Of course when the iterations are prematurely 
terminated, the locations of the photon absorptions 
and exits are not quite as accurate as when the itera- 
tions are continued, but the error is always less than 
the ratio of the flux that is remaining to that which 
has already escaped. In any case, using eq.(^|) is more 
intelligent than considering the remaining flux to be 
all escaped, all absorbed, or split up using just the 
scattering albedo. 

2.2. Radiative Transfer in Two-Phase Clumpy 
Media 

The main parameters defining a two-phase clumpy 
medium are the volume filling factor of the clumps, / c , 
the ratio of the clump to interclump medium (ICM) 
densities, a = p c / Picrm the total dust mass M, and 
the volume V of the medium. Typical ranges that 
may describe parts of the interstellar medium are 



0.01 < f c < 0.3 and 10 < a < 1 7 flSpitzer 1978|, |van 
Burcn 1989|, |Murthy et al. 1992|, [Gaustad fc van Bu 



ren 1993). The number of clumps and the clump sizes 
are secondary parameters, inversely related, which we 
shall discuss later. To obtain a formula for the ICM 
density, first define the density of the equivalent ho- 
mogeneous medium to be phom = M/V. Then since 
the total dust mass is simply the sum of mass in each 
phase, the homogeneous density can be expressed as 



Phc 



M c + Mi cm 
V 

VcMc Vicm M icm 

V V c V v lcm 

fcPc + (1 - fc)Pi cm 

{pc - Picm)fc + Picm , (7) 
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where M c and M icm are the dust masses in the clumps 
and ICM respectively, V c and Vi cm are the respective 
volumes, and so p c = M c /V c is the density of dust 



in the clumps, p^ 



M icm /V icm the density in the 



ICM, and f c = V c /V is the filling factor of the clumps. 
Equation (0) can be rearranged by substituting for p c 
with otpicm to give an equation for the ICM density 



Pho 



(a - 1) f c + 1 



(8) 



otphom from below, 
I a and p c — > p hom 



As f c — > 0, jOjcm — > p^orn and p c 
whereas if fc ^ 1, Picm -* Phom 
from above. 

One possible type of two-phase clumpy medium 
is that of cubic clumps randomly located on a body 
centered cubic lattice. The construction is similar to 
that used in percolation theory, where the probabil- 
ity of lattice site occupation is equivalent to the fill- 
ing factor and the clumps are then the occupied sites, 
and the cubic clumps do not overlap. Witt & Gordon 
(19 |96)| studied radiative transfer in such a two-phase 
clumpy medium defined on a cubic lattice, and we re- 
peated some of their Monte Carlo simulations giving 
a successful test of our methods. In this work we use 
spherical clumps defined on a high resolution 3D grid, 
allowing clump overlaps and the random locations of 
clumps to be chosen at a higher resolution than the 
clump size. In addition, the radiative transfer proper- 
ties are then readily approximated by analytic formu- 
las, as we demonstrate in §|| and §0. The stochastic 



media in the analytical m odels of |Boiss6 (1990)| and 
Hobson fc Scheuer (1993) do not impose any restric- 



tions on the shapes of the clumps, however, the mod- 
els are developed for the case of isotropic scattering. 

To construct a two-phase medium with randomly 
located clumps that may overlap, we must calculate 
the number of clumps needed to achieve the desired 
filling factor. This can be derived from the volume of 
a clump, v c , the total volume of the medium, V , and 
the filling factor, f c , but the possible occurrence of 
overlapping clumps complicates the calculation. As- 
sume that the clumps are all identical. Defining the 
volume fraction of a single clump as 



P = V 



(9) 



then the probability that a random point is not in any 
clump is equal to (1 — p) N °, were N c is the number of 
clumps, and therefore 



(see Appendix g for more discussion). Equation (10) 
is easily solved for the total number of clumps, 



Nr. 



Mi - /,) 

ln(l -p) 



(11) 



in terms of the total filling factor, f c , and the sin- 
gle clump filling factor, p. Although the above equa- 
tions apply to identical clumps of any shape, we 
shall consider only spherical clumps of radius r c with 
v c = 47rr^/3. Given f c ,r C7 and V, the clumpy medium 
is constructed by calculating N c using equations (||) 
and (fTl]), then selecting N c random points uniformly 
distributed in V, placing a sphere of radius r c and 
density p c around each point, setting the density be- 
tween the spheres to Pi C7n , and finally discretizing the 
medium on a high resolution 3D grid. When the ran- 
domly located clumps do overlap, the densities are 
not summed, so the medium is always characterized 
by only two possible densities: that of the ICM and 
the clumps. 

Monte Carlo simulations of radiative transfer were 
computed for three types of photon sources in a two- 
phase clumpy medium contained within a sphere of 
unit radius R$ — 1, in which the spherical clumps 
have radius r c = 0.05 and are 100 times denser than 
the ICM (a = 100). The 3D rectangular grid used 
to represent the medium has resolution of 127 3 vox- 
els, so that the clump diameters are about 6 grid 
elements. Figures [|, ^ and || show our results for 
an array of clump filling factors and optical depths, 
keeping r c and a constant. Three million photons 
were followed in each Monte Carlo simulation. Each 
of the images is the map of photons absorbed in a 
2D slice (127 x 127 pixels) through the center of the 
sphere. The grey scale from black to white indicates 
minimum to maximum absorption, on a logarithmic 
scale. Figure [j] represents the case of a single cen- 
tral isotropic luminosity source, analogous to a star 
or a centrally condensed cluster of stars in an H II re- 
gion. Figure || represents the case of isotropic photon 
emission distributed uniformly in the sphere, analo- 
gous to a uniform distribution of stars in a galaxy. 
Figure represents the case of an external uniformly 
illuminating source of isotropic photons, analogous to 
cold molecular clouds illuminated by the diffuse in- 
terstellar radiation field (IRF). In all cases the dust 
is characterized by a scattering albedo of uj = 0.6 
and asymmetry parameter g = 0.6, which are t ypical 



fc = l-(l-p) 



N,. 



ML 



values for UV p hotons scattering off dust grains ( Gor- 
don ct al. 1991 ). Recall that g = results in isotropic 
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horn 



= 5 



horn 



= 40 



Fig. 1. — Spatial distribution of photons absorbed 
by dust in a two-phase clumpy medium heated by a 
central point source. The simulations were performed 
by Monte Carlo calculations, and the images depict a 
2D slice through the center of the sphere. Scaling is 
logarithmic from the minimum of 0.1 photons (black) 
to the maximum of 4 x 10 4 photons absorbed per voxel 
(white). For more details see §2.2. 

scattering and g = 1 gives forward scattering. 

Moving vertically in the array of slices in the Fig- 
ures corresponds to increasing the volume filling fac- 
tor f c of the clumps while keeping the total dust 
mass constant (in each column). Thus the number 
of clumps appearing in each slice increases and pi cm 
decreases from bottom to top. Moving horizontally 
increases the equivalent homogeneous optical depth 
of the sphere, 



Tho 



Rs 



(12) 



which is the radial optical depth of extinction (absorp- 
tion plus scattering) that would result if the dust were 
distributed uniformly throughout the sphere instead 
of being clumpy. Increasing Thom can be viewed as 
either increasing the dust abundance (increasing pi cm 
and p c ) or changing the wavelength of the photons, 
either case resulting in more absorption. 

For the case of a central source (Figure |IJ), as 




horn 



t =A0 

horn 



Fig. 2. — Same as Figure [l] for dust heated by a 
uniform distribution of internal sources. Scaling is 
logarithmic from the minimum of 0.1 photons (black) 
to the maximum of 40 photons absorbed per voxel 
(white). For further description see §2.2. 
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r =A0 

horn 



Fig. 3. — Same as Figure [l] for dust heated by an ex- 
ternal uniformly illuminating source. Scaling is log- 
arithmic from the minimum of 0.1 photons (black) 
to the maximum of 80 photons absorbed per voxel 
(white). For further description see §2.2. 



Viom increases the ICM absorbs more photons and the 
clumps become opaque, creating the apparent shad- 
ows behind the clumps, when f c is small. However 
scattering by the dust causes photons to go behind 
the clumps and become absorbed, thus diminishing 
the effect of what would otherwise be completely dark 
shadows in the case of no scattering. As the clumps 
become opaque (-Thorn = 40), absorption occurs more 
at the clump surfaces. When f c is increased (keep- 
ing the total dust mass constant) the effect of shad- 
ows merges into the appearance of an absorption cav- 
ity. Similar effects occur for the case of uniformly 
distributed sources (Figure ||), however not as dra- 
matic as when the clumps are illuminated by a central 
source, since there are no shadows when the photons 
are emitted everywhere in the medium. At high f c 
the clumps dominate the medium and absorb most of 
the photons. In the case of a uniformly illuminating 
external source (Figure ||), when r^ om and f c are large 
we find that most of the impacting photons are ab- 
sorbed in the outer layers of the spherical region and 
the center is shielded from radiation. However, for 
small f c the center of the sphere is less shielded and 
instead the clump centers are shielded, most photons 
being absorbed on the surfaces of the clumps. In gen- 
eral, the clumpy medium allows photons to penetrate 
farther into the spherical region than if the medium 
was homogeneous, and this fact was the motivation 



and a conclusion in the work of Boisse (1990) , Hobson 



Scheuer (1993), and HP93, for the case of plane- 



parallel slab geometry. 

By observing the fraction of photons that escape 
in the case of a central point source we can compute 
the effective optical depth of the medium at a given 
wavelength A as 



T S (A) = - In 



-^esc(A) 



L (X) 



(13) 



where Lq(A) is the luminosity of the central source 
and Lesc(A) the escaping luminosity (for simplicity 
we shall hereafter drop any explicit dependence on A). 
We use the notation t$ instead of r e // because T e // 
is reserved for the case of no scattering (absorption 
plus scattering considered as interaction), whereas t$ 
includes the effects of scattering and is therefore de- 
pendent on the geometry of the medium, which in 
this case is spherical. Equation ( [l3"| ) is analogous to 
eq.([j]), and it follows that t$ < T e ff < Thom for a 
clumpy medium. Figure ^ shows the behavior of t$ 
as a function of Thorn for three cases of f c , with other 
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Fig. 4. — The effective optical radius, ts, of clumpy 
media in a sphere is shown as a function of the equiv- 
alent homogeneous optical depth, Thom, for three val- 
ues of f c , the clump filling factor. The symbols 
show ts given by eq.(|l^) using the results of Monte 
Carlo simulations including scattering (u> = 0.6 and 
g = 0.6), with f c = 0.1 (squares), 0.2 (triangles), and 
0.3 (diamonds). See §2.2 for more discussion. The 
curves are produced by analytic approximations de- 
veloped in this paper, and are discussed in §||. 



parameters having the values mentioned above. The 
diamonds indicate results of Monte Carlo simulations 
when f c = 0.3 (corresponding to top row in Figure [l]), 
the triangles show f c = 0.2, and the squares show 
f c = 0.1. There are two slopes in each curve of t$ 
versus r/j 0m , corresponding to the two phases of the 
medium. The steeper slope at low values of Thom is 
due to absorption by the increasing optical depth of 
the clumps. As Thom oo the clumps become opaque 
and any further increase in absorption is due to the 
ICM, and because it has a lower density the slope of 
the ts curves is less. The solid, dotted and dashed 
curves are produced by the analytic approximations 
presented in sections || and |[ and the agreement with 
Monte Carlo results is discussed further in «j. 



ESCAPE AND ABSORPTION 
HOMOGENEOUS MEDIA 



IN 



Since Monte Carlo simulations can require a large 
amount of computer time, it is useful to have ana- 
lytical approximations for the basic results of radia- 
tive transfer, such as, the fraction of photons that 
escape a bounded medium. In this section we present 
such approximations for homogeneous spherical me- 
dia with internal or external sources. The escape and 
absorption probability approximations are also com- 
pared with results of Monte Carlo simulations to as- 
sess their accuracy. These escape probability formu- 
lae are later applied to the case of clumpy media. 

3.1. Central Isotropic Point Source 

Consider an isotropic point source in the center of 
a spherical homogeneous medium. When the medium 
does not scatter photons the escape probability is sim- 
ply e~ T , where r = pnR is the optical radius of the 
sphere, characterized by a radius R and a mass den- 
sity p, and where k is the absorption cross-section of 
the dust per unit mass. If the medium also scatters 
photons, we can construct an analytical approxima- 
tion for the effective optical radius of the sphere, ts, 
defined as the negative natural logarithm of the escap- 
ing fraction of photons, as in eq. (|l3|) . The approxi- 
mation formula is based on limiting cases of the opti- 
cal parameters. When the optical depth of extinction 
(absorption and scattering), r = r a b s + T scat , is large 
and the scatt ering is isotropic (g ~ 0), then thcorct 



ical analysis flRybicky fc Lightman 1979, pp.37-38| ) 
suggests that ts is approximately 



t s ~ r (1 - uj)2 ; ( T > 1 and g ~ 0), 



(14) 



where u) = r SCQt /r is the scattering albedo. When the 
optical radius is small (t « 1) and g is any value, we 
expect that 



Ts 



r(l 



(t«1 and any g), 



(15) 



which is also an exact formula in the case of purely 
forward scattering at any optical depth. We can in- 
terpolate between these extreme cases using the fol- 
lowing formula: 

T S (T,u,g) = t(1-w)* t '»> , (16) 

where the interpolation exponent is given by 

X(r,9) = 1 - i(l-e-^ 2 )(l-# . (17) 
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Fig. 5. — Contours show the difference between the 
effective optical depth of a homogeneous sphere as 
predicted by eq.(|l6|) and that resulting from Monte 
Carlo simulations. 



The probability of escaping from the homogeneous 
sphere is then defined as 



T > t ac {T,w,g) = exp[-r 5 (r,w, ff )] 



(18) 



where the superscript "c" indicates that this is for a 
central point source. 

To check this approximation, the Monte Carlo code 
was used to compute a 3D matrix of effective optical 
radii as a function of optical parameters in the ranges 
< t < 13 with At < 1, < u < 1 with Alu = 0.1, 
and < g < 1 with Ag = 0.1, following more than 10 6 
photons in each simulation. Figure shows contours 
of the difference between ts(t, cj, g) given by eq.([i6|) 
and Ts from Monte Carlo simulations. The differ- 
ences are shown on 2D slices through the 3D parame- 
ter space, as a function of (a): (<?, lu) with r = 10; (b): 
(<7,t) with lu = 0.6; and (c): (w,r) with g = 0. The 
bold line in each panel depicts the region in parameter 
space where the agreement is perfect. The dotted con- 
tour lines indicate where eq.(|l6|) gives values less than 
the Monte Carlo results. The worst case of the ap- 
proximation occurs when uj = 0.6, g = 0, and r > 4, 
as panels (b) and (c) show that Ts(r,uj,g) underesti- 
mates the Monte Carlo results by a value exceeding 
unity. Also, panels (a) and (c) show that eq.(|l4|) is not 
always a good approximation. However the difference 
is relatively small compared to the value of ts « 10 
at (T,uj,g) = (13,0.6,0). In addition, when modeling 
radiative transfer through interstellar dust, isotropic 
scattering occurs at near-infrared (NIR) wavelengths 
for which the optical depth is smaller than at UV 
wavelengths where the scattering is more forward di- 
rected, and so the approximation is likely to be good 
at all wavelengths, as shown in panel (b). 

3.2. Uniformly Distributed Emission 

Consider a spherical homogeneous medium with a 
uniform distribution of isotropically emitting sources. 
When there is no scattering, Osterbrock (1989, pp. 385- 
386) derived an exact solution for the photon escape 
probability given by 



1 



1 

2t2 



1 

272 



(19) 



where r = pnR is the optical radius of the sphere 
(R = radius of sphere, p = density of absorbers, and 
k = absorption cross-section). We find that P e {r) 
agrees exactly with the escaping fraction of photons 
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computed in Monte Carlo simulations. Using the op- 
tical depth of extinction resulting from both absorp- 
tion and scattering, r ext = r abs + T scat , we can get 
the probability, P e {T ext ), that photons will escape the 
spherical medium without any interactions (absorp- 
tion or scattering), and we call it the extinction es- 
cape probability. Making the assumption that the 
scattered photons have the same spatial and angular 
distribution as the emission sources, the extinction es- 
cape probability can be applied recursively to arrive 



at the scattering escape probability |(Lucy ct al. 1991 



Pe(r) 



1-«[1-P e (r)] 



(20) 



where r = T ext and oj — T scat /T ext is the scattering 
albedo. This formula estimates the effects of scatter- 
ing, namely that scattered photons may eventually es- 
cape, thereby increasing the escape probability. The 
superscript u u" indicates that "P" sc is for uniformly 
distributed emitters in a sphere. Actually, Pe(r) in 
the formula can be any extinction escape probability 
for any geometry, but this paper will concentrate on 
the case of spherical geometry, for which eq.(19) ap- 
plies. The derivations of equations (19) and (2C) are 
reviewed in Appendix 

Equations (19) and (pp|), which we shall call the 
Osterbrock-Lucy escape probability (OLEP) formula, 
were tested extensively against Monte Carlo radia- 
tive transfer simulations with multiple scattering and 
found to be a reasonable approximation of the fraction 
of photons escaping from a spherical homogeneous 
medium. Since the scattering asymmetry parameter, 
g = (cos 9 Bca t}, does not enter into the formula for 
V^ sc {t, uj), the value of g for which the OLEP formula 
agrees exactly with Monte Carlo simulations is found 
to depend on the extinction optical depth r. The up- 
per panel in Figure ^ shows contours of the percent 
relative difference between the OLEP formula and es- 
caping fraction of photons found by Monte Carlo cal- 
culations (following 10 6 photons in each simulation), 
as a function of g and r, when the albedo is uj = 0.7. 
The values of g for which the OLEP agrees with 
Monte Carlo results are indicated by the zero contour 
level drawn with a thick solid line. The dotted con- 
tours indicate where the escape probability formula 
underestimates the escaping fraction, and the thin 
solid contours indicate where it overestimates the es- 
caping fraction. For optically thin situations (r < 1) 
the OLEP formula agrees well with the Monte Carlo 
simulations for all values of g, with the best agree- 



h 10 




i- 10- 



Fig. 6. — Analysis of the escaping fraction of photons 
that are emitted by sources uniformly distributed in 
a homogeneous sphere, as a function of the scattering 
asymmetry parameter g, and r, the extinction optical 
radius of the sphere, with albedo uj = 0.7. The top 
panel shows contours of the percent relative difference 
between the escaping fraction given by the OLEP for- 
mula [eqs.(|l9|) and (pp|)1 and that predicted by Monte 
Carlo simulations. The bottom panel shows contours 
labeled with the escaping fraction predicted by Monte 
Carlo simulations, over-plotted with the zero differ- 
ence contour (thick solid line) from the top graph, 
referred to as g*(j). The dashed line is an approxi- 
mation of g* (t) given by eq. (^ll) . 
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merit occurring when the scattering is nearly isotropic 
(g ~ 0) . As the optical depth increases the agreement 
shifts toward more forward scattering cases (g — > 1) of 
the Monte Carlo simulations. The maximum relative 
difference in the range shown is 26% overestimation 
at t = 20, lu — 0.7, and g = 0, where the escaping 
fraction is only 0.10. The behavior of the relative dif- 
ference is similar for other values of the albedo, and 
always decreases with albedo, to zero as lu — > 0, since 
then V^ sc (t,lu) — > P e {T~). The value of g for which 
V^ sc (t,lu) agrees with Monte Carlo calculations (the 
zero difference contour) is found to be a function of 
optical depth only, independent of the albedo lu, and 
so we designate this special value by g*(r). 

The lower panel in Figure ^| shows contours la- 
beled with the escaping fraction resulting from the 
Monte Carlo calculations, and the thick solid line is 
again g*(r), the zero difference contour from the up- 
per panel. We find that a reasonable empirical ap- 
proximation for g*(r) is given by 



purely forward scattering: 



9*{t) 



1 - 



3.3 



3.3 



(21) 



as shown by the dashed line. 

There is a simple explanation for why the best 
agreement between P™ sc (t,lu) and Monte Carlo sim- 
ulations occurs at a value of g that approaches unity 
as t — ► oo. To explain this fact we rewrite eq. (p(i|) to 
read 



(22) 



lu + (l-w)/P e (r) ' 
where P e (T~) is the extinction escape probability from 



a sphere given by eq.(19). As r — > oo the extinction 
escape probability approaches zero as P e (r) ~ 3/4r . 



(23) 



Applying this limiting behavior to eq.(|22[) gives 



lu + 4t(1-lu)/3 4t(1-w) 



since the term 4t (1 — lu)/3 3> uasn oo. Now we 
recognize that r (1 — u>) is exactly the optical depth 
in the case of purely forward scattering, when 5 = 1, 
since then scattering acts like a reduction in absorp- 
tion. The forward scattering escape probability is 
then exactly P e [r(l — lu)] , which has the same limit- 
ing behavior as eq.(B3J) when r — * oo : 



Pe [T(l 



4t(1-u;) 



(24) 



Therefore as r — * oo the OLEP approximation ap- 
proaches from below the exact escape probability for 



V: sc (r^) - PeHl-Lu)} 



(25) 



Since the left side is valid at g — g*(r) and the right 
side is valid at g = 1, the limiting behavior implies 
that g*(r) — > 1 as r — > oo, which is satisfied by the 
approximation equation (^l|). 

In summary, V^ sc (t, lu) overestimates the actual es- 
cape probability when g < g*(r) and underestimates 
it when g > <?*(t). For values of g in the range 
g*(r) < g < 1 the escape probability approximation 
can be improved by linear interpolation in g between 
the values V^ sc (t,lu) at g*(r) and P e [t(1 — lu)] at 
g = 1. If an analytical approximation for the escape 
probability when g = — 1 (purely backward scatter- 
ing) were known, then three point parabolic interpo- 
lation could be used to obtain a very good approxi- 
mation to the actual escape probability for all g, and 
hence all (r,u>,g). When modeling absorption and 
scattering by dust, large optical depths usually occur 
at UV wavelengths for which the scattering by dust 
tends to be more forward oriented, having values of g 
near g*(r), so the error incurred by using the OLEP 
approximation may be relatively small. 

3.3. Finite Number of Uniformly Distributed 
Point Sources 

The question arises as to whether the escape prob- 
ability formula can be successfully applied to a dis- 
crete collection of randomly placed isotropic point 
sources ( "stars" ) instead of a continuous uniform dis- 
tribution of emitters, and in particular what is the 
minimum number of stars for which the escape prob- 
ability formula gives correct results. To answer these 
questions we performed Monte Carlo simulations for 
cases of 2 to 80 randomly placed "stars" in a homo- 
geneous sphere, and for each case we simulated 200 
trials, where in each trial 10 5 photons were followed to 
find the fraction that eventually escapes. The results 
are shown in Figure for when the optical radius of 
the sphere is r = 5, for the case (a) of no scatter- 
ing (lu = 0) and case (b) with scattering (lu — 0.6 
and g — 0.6). The horizontal lines in cases (a) and 
(b) are the escape probabilities P e (r) and T 7 " e (r, m), 
respectively, given by equations ( |l9| ) and (pGJ). Each 
diamond symbol shows the average escaping fraction 
of the Monte Carlo trials for a fixed number of stars, 
and the error bars are plus and minus one standard 
deviation from the average value. The distribution 
of escaping fractions from the Monte Carlo trials is 
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Fig. 7. — Monte Carlo simulations of the fraction of 
photons that escape from a homogeneous sphere of 
dust after being emitted by a small number of point 
sources ("stars"), as indicated on the horizontal axes. 
Case (a) is for absorption only (u = 0), case (b) is 
with scattering (lj = 0.6 and g = 0.6), and the opti- 
cal radius of the sphere is r = 5 in both cases. The 
diamonds and error bars are each the means and stan- 
dard deviations of 200 Monte Carlo simulations. The 
solid horizontal lines are predictions by the OLEP 



formula. The figure is discussed in §3.3 



approximately Gaussian. We find that the standard 
deviation of the trials goes as the inverse square root 
of the number of stars n s , as expected for sums of 
random variables, and in particular <j ~ 0.1 (n s )~ - 5 
for all r > 1 . Since the escaping fraction in the case of 
no scattering behaves as P e {j) ~ 3/4r when r — > oo, 
the expected relative deviation (a/P e ) in the escaping 
fraction of photons from n s stars randomly located 
in a sphere of optical depth r > 1 will be approxi- 
mately 0.133 ttt-jT 0,5 for the case of absorption only, 
and less when scattering is also involved. Thus the 
number of stars must increase as the square of the 
optical depth to maintain the same expected relative 
deviations. The analytic escape probability agrees 
well with the average escaping fraction found for each 
group of Monte Carlo trials, it is just a question of 
what relative deviation is acceptable. 

3.4. Uniformly Illuminating External Source 

In this section we introduce approximations for the 
probability that externally emitted isotropic photons 
will interact with and get absorbed by a homogeneous 
medium contained in a sphere. The term interaction 
includes both absorption and scattering events. Given 
that a distribution of photons encounters the sphere 
at a random impact parameters, the probability of 
interaction at any location inside the sphere is 



1 



2r z \ r 2r z 



1 



(26) 



where t = pnR is the extinction optical radius of 
the sphere. This equation is derived by averaging the 
transmission over all possible impact parameters and 
computing the ratio of net absorption versus impact- 
ing flux, using the same techniques as for the deriva- 
tion of eq. (|l9|) (see Appendix |c]). When t < 1 we 
have Pi{r) ~ 4r/3, which is the optical path length 
through a sphere averaged over all i mpact parame- 



ters. Equation (^6|) was utilize d by Neufcld (1991 



and Hobson fc Padman (1993) in their models for 



clumpy media, as we shall discuss in §|llj. Note that 
equations (26) and ( |l9| ) are related by 



A(T) = yP e (T) 



(27) 



giving a duality between absorption (interaction) of 
an external source and escape (non-interaction) of an 
internal source. 

If the medium scatters photons, we would like to 
know what fraction of Pi(r), the interacting frac- 
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tion, will eventually escape due to multiple scatter- 
ing events, and thus determine what fraction of pho- 
tons are actually absorbed. In the special case of 
purely forward scattering there is an exact formula 
for the absorbed fraction of photons. When the scat- 
tering deflection angle is always zero, it is equivalent 
to no scattering with a reduced optical depth equal 
to t(1 — w), where oj is the scattering albedo. So 
the actual fraction of photons absorbed in the purely 
forward scattering case is 



(28) 



where the superscript u x" indicates that this is for 
the case of an external source and g is the scattering 
asymmetry parameter. As expected, equations (^) 
and ((H) agree exactly with Monte Carlo simulations. 

For non- forward scattering cases, when g < 1, we 
can use the methods discussed in S3. 2 for the case of 



internal uniformly distributed emitters, namely the 
OLEP equations (|l9| ) and (|2Cj). Given that photons 
interact with the medium inside the sphere, the prob- 
ability that they scatter is simply ui. These scatter- 
ing events can be regarded as re-emitted photons, and 
assuming that they are approximately uniformly dis- 
tributed in the sphere, the scattered photons have a 
probability V" sc (t,uj) of escaping [eq.(|20|)]. The frac- 
tion of the interacting photons that actually get ab- 
sorbed is then 1— luV" sc (t,u>). Thus an approxima- 
tion for the fraction of photons absorbed is 



V x abs {r,u) = Pi(r) [1-u>VZJt,u>)\. 



(29) 



Recall that V^ sc (t,ui) is valid only for g = g*(r) as 
discussed in § |3.2| , and so we expect "P„ bs (r, oj ) to fol- 
low the same pattern. When r c is large, the pene- 
tration depth of photons decreases so they will not 
be uniformly distributed in the sphere, as assumed 
above, but this does not severely affect the applica- 
tion the approximation as we show next. 

The upper panel in Figure || shows contours of the 
percent relative difference between V^ bs (T,uj) and ab- 
sorbed fraction of photons found by Monte Carlo cal- 
culations, as a function of g and t, for lo = 0.7. The 
values of g for which eq. (^) agrees with Monte Carlo 
results are indicated by the zero contour level drawn 
with a thick solid line. The dotted contours indicate 
where the absorption probability formula ( p9| ) under- 
estimates the absorbed fraction, and the thin solid 
contours indicate where it overestimates. The dif- 
ference is small and independent of g for optically 
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Fig. 8. — Analysis of the absorbed fraction of pho- 
tons from an external source which uniformly illu- 
minates a homogeneous sphere, when the scattering 
albedo is u> = 0.7. The top panel shows the percent 
relative difference between the absorption probability 
given by eq.(p9|) and the absorbed fraction found by 
Monte Carlo simulations, as a function of the asym- 
metry parameter g, and r, the extinction optical ra- 
dius of the sphere. The contours in the bottom panel 
are labeled with the absorbed fractions computed by 
Monte Carlo simulations. The thick solid line is the 
locus of agreement between V% bs [eq.(^9])] and Monte 
Carlo simulations (zero difference contour from top 
graph), referred to as <?*(t). The dashed line is an 
approximation of g*(r) given by eq.(|2l|). 
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thin situations (r < 1), and as the optical depth 
increases the agreement shifts toward more forward 
scattering cases (g — > 1) of the Monte Carlo simula- 
tions. Not surprisingly, the behavior of the relative 
difference contours is identical to that in Figure o 



of §3.2, since we are using V^ sc {t,uj) in the formula. 
The relative difference decreases as ui — ► 0, since then 
r Pabs( T ' u> ) ~~ * ft ( T ) an d Pi ( T ) i s an exa ct formula. 

The lower panel in Figure |^ shows contours labeled 
with values of the absorbed fraction obtained from 
Monte Carlo simulations of a uniformly illuminating 
isotropic external source. Comparison with Figure ^ 
shows that the behavior of the absorbed fraction of 
an external source is quite different from the absorbed 
fraction of an internal source (one minus the escaping 
fraction shown in lower panel of Figure ^) . Note that 
for g < 1, as r — > oo the absorbed fraction tends 
asymptotically toward a constant value, since in the 
opaque limit scattering occurs at the surface of the 
sphere. The thick solid line is again the contour of 
zero difference between V^ bs (r, ui) and Monte Carlo 
simulations, and the dashed line is the approximation 
of g*(j) given by eq. (|2l|) . This zero difference contour 



is found to be independent of u> and fit well by eq.(21), 
also because V^ sc (t,uj) is used in the approximation 
formula. 

Let us now examine the behavior of the absorption 
probability as r — > oo to show in more detail why the 
best agreement with Monte Carlo simulations occurs 
at a value of g that approaches unity. Using eqs.(|27|) 
and (ptj| ), and defining P = P e (r) for convenience, 
eq.(|29|) can be rewritten as 



3 

lip 

3 

4r(l 



1 - 



l-w(l-p) 

1 - U 



1-W(1-P) 



P 



3 

4r(l-co) 



l-w(l-P) 



(30) 



thus extending eq.(p7|), the duality between the cases 
of internal and external sources, to the case of lu > 0. 
Equation (p8f) can be rewritten as 



V* bs (T,Lu;g = l) = 4r(1 w) P e [T(l-u,)] (31) 



where we have used eq.(|27|) again. The analysis given 
by eqs.(HH|) proved that T? sc {t, lu) -> P e [r(l - u)] 



from below as r — > oo, thereby proving that 

r: bs (r, lu; 9 = g*) -> V x abs (r, u;g=l) (32) 
from below, implying that g*(r) — > 1 as r — > oo. 



4. ESCAPE AND ABSORPTION 
CLUMPY MEDIA 



IN 



A viable approach to the problem of estimating 
the escaping and absorbed fractions of photons in 
an inhomogeneous medium is to find effective values 
for the absorption and scattering properties of the 
inhomogeneous medium and then use these effective 
values in the escape/ absorption probability formulae 
that were developed for homogeneous media. The 



work of Hobson & Padman (1993) (HP93) provides 
the means to such an approach for two-phase clumpy 
media by a model they called the "mega-grains" ap- 
proximation, which reduces the inhomogeneous radia- 
tive transfer problem to an effectively homogeneous 
one. In the following sections we review the mega- 
grains approximation, and then discuss how it is uti- 
lized to give escape/absorption probability approxi- 
mations for clumpy media. We also introduce some 
improvements to the mega-grains approximation of 
HP93. In addition, we present new formulae giving 
the approximate fraction of photons absorbed in each 
phase of the medium, which can then be used to es- 
timate the dust temperature and infrared emission 
from the clumps and ICM separately. 

4.1. The Mega-Grains Approximation 

Natta & Panagia (1984) developed analytic ap- 
proximations for the extinction by dust having var- 
ious kinds of inhomogeneous distributions, in par- 
ticular, a clumpy distr ibution with an empty inter- 
clump medium (ICM). Neufeld (1991) proposed the 



approach of treating spherical clumps in a two-phase 
medium as large grains, and applied it to the special 
case of the scattering, absorption, and escape of Lya 
radiation, involving both gas and dust. HP93 pro- 
posed a more general set of formulae called the mega- 
grains approximation, which also treats the spherical 
clumps as large grains but with absorption and scat- 
tering coefficients in direct analogy with dust grains. 
Their mega-grains model gives equations for the ef- 
fective optical depth and effective scattering albedo 
of a two-phase clumpy medium, with a non-empty 
ICM. In this section we review the derivation of 
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the mega-grains approximation and introduce an im- 
proved equation for the effective scattering albedo. 
We also develop a new approximate formula for the ef- 
fective scattering asymmetry parameter of the clumpy 
medium. To distinguish between similar properties of 
the clumps (mega-grains) and the dust we will assign 
them the subscripts "c" and "d", respectively. So lu c 
will be the macroscopic albedo of the clumps and iOd 
is the microscopic albedo of the dust. Also, an im- 
plicit dependence on wavelength is assumed in all the 
following. 

4.1.1. Effective Optical Depth 

Assume that a collection of spherical clumps, each 
of radius r c and having dust mass density p c , are ran- 
domly distributed in an inter-clump medium (ICM) 
ha ving a constant dust mass density of p^m < o r . 



Fol lowing HP93. define the clump optical radius 



\Pc Picm) K r c 

[a - 1) np icm r c 



(33) 



where n is the total absorption plus scattering cross- 
section of the dust grains per unit mass and recall that 
a = p c j picm- The reason for the subtraction of densi- 
ties in this definition is to disregard the clumps in the 
limit of when the clump density approaches the den- 
sity of dust in the ICM, since then there are effectively 
no clumps. Thus the mega-grains are considered to be 
just the over-density, p c —Picm, allowing for analysis as 
a separate component superimposed on a continuous 
medium of lower density. Now considering the clumps 
as large grains, the probability that a randomly emit- 
ted photon will encounter a clump is just irr^.. Once 
a clump is encountered, the probability that the pho- 
ton will interact (get absorbed or scatter) with dust 
in clump is Pi(r c ), given by eq.(|26|), the interaction 
probability for a sphere of optical radius r c . Then we 
can define the interaction coefficient per unit length, 
A mg , in the medium for just the mega-grains as 



A, 



n c -Kr 2 c Pi{r c ) 



(34) 



where n c is the number of clumps per unit volume. 
Here the use of the interaction probability Pi(r c ) is 
analogous to the absorption and scattering efficiency 
coefficient of dust grains. 

Define the effective interaction coefficient per unit 
length, A. e ff, in the two-phase clumpy medium by 
including the ICM: 



A 



eff 



Kpi, 



(35) 



This reduces to the interaction coefficient for a ho- 
mogeneous medium, Kphom, when a — p c / Picm = 1 
or f c = 0, i.e. when the medium is homogeneous 
[see eqs.(fj]) & (f|)]. For a clumpy medium in a plane- 
parallel slab of thickness L, or sphere of radius R = L, 
the effective optical depth is then 



T eff 



LA 



eff 



LA r 



Lnpi, 



(36) 



where by effective we mean the optical depth corre- 
sponding to the average transmission of a large col- 
lection of randomly chosen paths through the slab, as 
expressed by eq.rtjj). In §4.2 we analyze in more de- 



tail the implementation of the mega-grains equation 
( |34| ) , developing an improved version, and discuss the 
behavior of r e // as a function of all the parameters 
characterizing the clumpy medium. Equation (|34|) is 
compared with the ap prox imation developed by |Natta 
(1984)| in §|J. 



Panagk 



4-1.2. Effective Scattering Albedo 

The effective scattering albedo of the two-phase 
clumpy medium is logically defined as a weighted 
combination of the albedos for each phase in the 
medium: 



A 



J eff 



my 



KPicm 



A 



eff 



(37) 



where u>d is the albedo of the dust and u> c is the effec- 
tive albedo of a clump, which we discuss below. First 
note that in the limit as the medium becomes homoge- 
neous (f c — > 0) we have A mg — > and A e // — * Kpi cm , 
so then u) e ff ~* ^>d as expected. 

Given that interactions with dust in a clump have 
occurred for a group of photons, the effective albedo of 
the clump is the fraction of those photons that man- 
age to eventually escape from the clump by means of 
multiple scattering. HP93 suggested an approxima- 
tion for the clump albedo for which we shall use the 
symbol uj^ p : 



.hp _ 



1 + (l-w d )4r c /3 



(38) 



where r c is the optical radius of a clump as given by 
eq.(|3|). As r c — * we have uif/~ p — * u>d as required. 
As t c — > 00 then oj^ p — > 0, but this cannot be gener- 
ally true since we expect always some backscattering 
from the surface layer of the clump as long as the 
dust scattering asymmetry parameter gd < 1. Note 
that the approximation does not consider the distri- 
bution of scattering angles. In addition we shall show 
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that for small r c , equation ( p8[ ) is below the minimum 
possible value for the clump albedo. 

Let us examine the special case of purely forward 
scattering, when the asymmetry parameter of the 
dust is 9^ = 1 . In this case there is an exact formula 
for the clump albedo. The forward scattering case for 
photons randomly impacting a spherical clump was 
discussed previously in the presentation of equations 
(26) and (p8|), and the actual fraction of photons ab- 
sorbed is Pi [t c (1 — ujd)] ■ The fraction that escapes is 
then Pi(T c ) — Pi [t c (1 — ujd)] , and the ratio of escaping 
fraction over interacting fraction is the clump albedo: 

" c = WI) (39) 

where the superscript "1" indicates that ga = 1- As 
t c — ► 0, we have Pi(r c ) — > 4r c /3 and then 

x _ 4r c /3-4r c (l-o; d )/3 
" c ^ 4r c /3 

= 1 - (1 - Vd) = OJd 

reaching the correct limit. On the other hand, if t c — > 
oo then u\ — > 0, which is expected in this case of 
purely forward scattering. We shall show that uj\ is 
the minimum clump albedo over all — 1 < gd < 1. 

To determine the actual clump albedo and how it 
depends on dust scattering parameters we performed 
Monte Carlo simulations of photons randomly im- 
pacting a spherical clump for many values of the 
optical parameters in the ranges of < t c < 25, 
< u>d < 1) and <«;<£< 1. The fraction of pho- 
tons that interact with the clump is found to be ex- 
actly given by eq.(p6|), as expected. Of these inter- 
acting photons, the fraction that scatters and even- 
tually escapes from the clump is the clump albedo 
uj c . Figure ^ shows u> c versus r c for three values of 
the dust albedo uid = 0.3, 0.6, and 0.8 (apparent 
when t c = 0). The square symbols indicate the case 
of gd = (isotropic scattering), the diamonds show 
gd = 1 (forward scattering), and intermediate values 
of < gd < 1 at increments of 0.1 are plotted as dots 
vertically connecting the squares and diamonds. The 
dotted line is ui\ given by eq.([39|), the dashed line is 
given by eq.(J38|), and the solid line we shall in- 



troduce shortly. Clearly, the theoretical forward scat- 
tering clump albedo, lo\ (dotted line), agrees exactly 
with the Monte Carlo results for <jy = 1 (diamond 
symbols), as expected, and is a lower bound on the 
clump albedo. From the dashed line it is apparent 



HP 



3 0.4 




Fig. 9. — Comparison of theoretical approximations 
for the clump albedo given by HP93, [ui^ p of eq. (|38|) , 
dashed line], and those developed in this paper, to* 
[eq. (|40|) , solid line] and uj\ [eq. (|39"|) , dotted line], with 
Monte Carlo simulations (squares for gd — 0, dia- 
monds for gd = 1), as a function of the clump optical 
depth r c . 
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that ujH p < u)\ when r c < 2, which violates the lower 
bound. 

Another formula for the clump albedo can be de- 
rived using the OLEP formula [eqs. ( |l9| ) and (pp|)] in 
the same manner as was discussed in the derivation of 
equation (|2"s|). Given that a photon interacts with the 
dust inside a clump, the probability that it scatters is 
just ujd- These scattering events can be regarded as 
re-emitted photons, and assuming that they are ap- 
proximately uniformly distributed in the clump, they 
have a probability P" sc (r c , w^) of escaping, thus ob- 
taining the clump albedo 



OJ* = UJ d V^ c (T c ,UJ d ) 



(40) 



which gives the solid line shown in Figure §. The 
escape probability equation ( p0| ) for V^ sc (T c ,uj d ) can 
be rearranged to give 



UJ,.. 



oj d 



W d + (1 - UJ d ) I Pe{r c ) 



(41) 



where P e (T c ) is the extinction escape probability, as 
given by eq. (p~9|) . As r c — > 0, we have P e (r c ) — > 1 and 
so uj* — ► idd, because after a single scattering a photon 
will most likely escape an optically thin clump. As 
t c — ► oo , we have P e (r c ) — > 3/4r c and then uj* — > 0. 
However, uj* gives the actual clump albedo for a single 
value of g d for each r c , and this value, g* d — g*(r c ), 



is approximately given by eq.(21), since we are using 
the approximate escape probability 'P" sc (t c , o^) (see 
Figure ^ . Therefore the fact that uj* — > as r c — > oo 
is just a consequence of the fact that g*(r c ) — > 1 as 
t c — > oo, coupled with the fact that the exact equation 



(39) for the case of g d = 1 gives the limit of zero for 
the albedo of increasingly opaque clumps. 

The clump albedo depends more strongly on gd as 
t c increases (see Figure ||) , creating a large spread in 
the values of uj c between the cases of forward scat- 
tering (gd — 1, diamonds) and isotropic scattering 
(gd = 0, squares). This is because as r c increases 
the average photon penetration depth to first scatter- 
ing decreases, and then forward scattered photons are 
likely to be absorbed in the clump whereas the spher- 
ical geometry presents many opportunities for the es- 
cape of any non- forward scattered photon.. Thus a 
higher probability of backscattering, corresponding to 
smaller values of gd, increases the probability that 
a photon will escape, giving a larger effective clump 
albedo. For any gd < 1 our Monte Carlo simula- 
tions indicate that the albedo of increasingly opaque 
clumps (t c > 20) remains non-zero. 



In summary, Figure ^ compares the theoretical ap- 
proximations for the clump albedo developed in this 
paper, uj* [solid line, eq.([H|)], uj\ [dotted line, forward 
scattering, eq.(|3^)], with the approximation given by 
HP93, wf p [dashed line, eq. (§§)], and with the Monte 
Carlo results (symbols). From the figure it is clear 
that uj* > uj^ p and uj* > uj\ for all clump optical 
depths, and for r c < 2 our proposed eq.(f4l|) for uj* 
docs better at predicting the clump albedo computed 
by Monte Carlo than uj^ p given by eq.(|38|). We shall 
assume that uj* given by eq.(^) is a sufficiently good 
approximation of the actual clump albedo, and will 
use it in eq. ( |37|) to obtain the effective albedo of the 
clumpy medium. If g*(r c ) < gd < 1, one can inter- 
polate in the variable gd between values uj* and uj\ 
corresponding to the endpoints of the range to obtain 
a more accurate estimate of the clump albedo. 

4-1-3. Effective Scattering Asymmetry Parameter 

Hobson and Padman did not consider the angular 
scattering distribution of photons in their paper in- 
troducing the mega-grains approximation, but a sim- 
ilar approach as used for the effective albedo of the 
clumpy medium can be followed. We define the ef- 
fective phase function asymmetry parameter for the 
two-phase clumpy medium as a weighted combina- 
tion of the asymmetry parameters of each phase in 
the medium: 



Jeff 



g c A r 



gd ^Picm 



A 



eff 



(42) 



where gd is the asymmetry parameter of the dust 
and g c is the clump scattering asymmetry parame- 
ter, which we discuss below. Note that g e ff — > gd 
as the medium becomes homogeneous (f c — * 0), since 
then A mg — > and A eff — * np icm . 

Let us assume that a collection of photons enter a 
clump with parallel ray paths, are scattered by the 
dust in the clump, and eventually escape from the 
clump. How does the exiting angular distribution de- 
pend on the dust scattering properties and the optical 
depth of the clump? The results of our Monte Carlo 
simulations shown in Figure [lO] answer this question. 
We computed g c — (cos 8 exit) for all simulations, 
where 9 ex u is the exiting angle respect to the enter- 
ing parallel rays, and plot g c versus the clump optical 
radius r c , for values of gd — {—0.4, 0.0, 0.6, 0.9} (ev- 
ident when t c — 0), and albedos 0.1 < UJd < 0.9. 
The square symbols are for when UJd — 0.1, the dia- 
monds for ujd = 0.9, and intermediate values of UJd in 
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increments of 0.1 are plotted as dots connecting the 
squares and diamonds. The distribution of exit an- 
gles does not actually follow the HG phase function 
with g — g c , since the scattering is complicated by the 
spherical geometry, however, the HG phase function 
is a reasonable approximation. In the case of forward 
scattering (gd > 0), the clump scattering distribu- 
tion becomes more isotropic as the optical radius r c 
of the clump increases and as the dust albedo u>d in- 
creases. If gd = (isotropic scattering by dust), then 
as the clump becomes opaque the clump asymmetry 
parameter approaches a backscattering distribution of 
g c = — 1/3 independent of the dust albedo. In addi- 
tion, when \gd\ is near zero and r c is large, it is ap- 
parent that g c approaches —1/3 as idd increases. Our 
computations of g c agree with the three cases com- 



puted and mentioned by Code fc Whitney (1995) 



The solid lines in Figure 10 are produced by the 
following empirical formula: 



9c(T c ,uJd,gd) = 9d - C 1 - 



1 



,-B/A 



1 + e (j.-B)/A 



(43) 



where 



A = 
B = 

C EE 



1.5 + 4g d + 2uj dy /g2exp(-5g d ) 

2 - g d (l - g d ) - 2u> d g d 

1 

3 - \f2~g~d - 2Lu d gd{l - gd) 



The formula is a good approximation of the Monte 
Carlo results when gd > 0. For gd slightly negative 
one can shift the gd = curve downward to get an ap- 
proximation of g Cl or just use g c = gd when gd < —0.2. 
The exact value of g c is of importance mainly for the 
case of a central point source, since the escaping frac- 
tion is then very sensitive to the distribution of angu- 
lar scattering. 

4.2. The Extended Mega- Grains Approxima- 
tion 

The mega-grains approximation as presented by 
Hobson & Padman is limited to small values of the 
clump filling factor. In addition, we found that for 
optically thin situations it predicts an effective opti- 
cal depth that is slightly greater than the equivalent 
homogeneous optical depth, in violation of eq.(^|). In 
this section we develop an improved version of the 
mega-grains approximation that resolves these prob- 
lems, and extends the approximation to all values of 
the clump filling factor. 




1 



1 5 



2G 



Fig. 10. — The effective clump scattering asymmetry 
parameter g c = (cos 9 exit) ^ where 9 e xit is the angle be- 
tween impacting and exiting directions as simulated 
by Monte Carlo calculations, is plotted as a function 
of the clump optical depth r c , for four values of the 
dust scattering asymmetry parameter, gd = —0.4, 0.0, 
0.6, and 0.9, corresponding to the intersection of each 
curve with the vertical axis at r c = 0. The square 
symbols are for dust scattering albedo of uj d = 0.1, 
the diamonds for ujd = 0.9, and the dots show in- 
termediate values of u> d in increments of 0.1. The 
solid lines show the empirical approximation given by 
eq.@ for uj d = 0.5 . 
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An important term in eq.(34) is n c , the density 
of clumps, which we shall now discuss in detail since 
it is the key to our extension of the mega-grains ap- 
proximation. Define the porosity, Q c , of a randomly 
located collection of identical clumps as the ratio of 
the total volume of clumps, including possible over- 
laps, to the volume V of the medium: 



V 



(44) 



where N c is the total number of clumps and v c is the 
volume of just one clump. Equation ([II]) is easily 
solved for the density of clumps 



N c Q c 3Q C 



V v c 4-7rrJ? 



(45) 



as a function of the porosity and the radius of a spher- 
ical clump r c . Substituting for n c in eq.(|34|) gives 



A mo — 



3Q C 
4rv 



Pi(Tc) 



(46) 



for the interaction coefficient of the mega-grains, 
where r c is the optical radius of a clump. 

The filling factor f c is related to the porosity, Q c , 
of the clumps by 



fc = 1 - e 



(47) 



because the probability that a random point is not 
in a clump goes as e~ < ^ c (this is a very good approx- 
imation to eq.([l0|) when v c <C V, see Appendix [^). 
When f c <C 1 then f c ~ Q c , however as f c — > 1 we 
have Q c — — ln(l — f c ) — > oo, since then the clumps 
tend to overlap. Thus an obvious problem with equa- 
tion (|^) is that as f c — * 1 we have Q c — > oo causing 
A mg — » oo. Indeed, Hobson & Padman found that 
/ c < 0.3 is the useful range for the mega-grains ap- 
proximation. The clump overlapping fraction is calcu- 
lated as (Q c — f c )/Qc, and for f c < 0.3 it is less than 
16% of the volume of clumps, but as f c > 0.3 the 
overlapping volume fraction increases rapidly toward 
100% so that the clumps can no longer be treated as 
separate mega-grains. 

There is another problem with the previous equa- 
tion, which is more subtle but still important. Con- 
sider the behavior of eq.(^) as t c — > 0, which occurs 
when either r c — > 0, k — ► or phom —> (hold- 
ing f c and a — p c / picm constant). In that case 
Pi(i~ c ) ~ 4t c /3, which simplifies eq.(Eq), and upon 



substituting for r c with eq.(|33|), the clump radius can- 
cels giving 



A ~ ^ c _ I 
l*-mg ~ 7~ c \Pc 



i) Qc 



(48) 



when the clump optical radii are very small. Substi- 



tuting this into eq.p5[) gives 

Ae// 



A 



mg T" fcpicm 
K . (/•*'- Picm)Qc 



Pi, 



(49) 



for the behavior of the effective interaction coefficient 
of the two-phase clumpy medium as r c — * 0. The 
problem with this behavior is that since Q c > fc we 
get that 



A e ff > K [(pc - Picm)fc + Pu 



Kpho 



(50) 



using eq.(0) for the final equality. Thus as the clump 
optical depths become small, A e // exceeds the in- 
teraction coefficient of the equivalent homogeneous 
medium. If L is the geometrical thickness of the 
medium then eq.(J50|) implies that 



" LA-eff > Lnp hom — T h , 



(51) 



contradicting equations (S) and (^) , which state that a 
clumpy medium is more transparent than the equiva- 
lent homogeneous medium. The inequality in eqs.(^Oj) 
and (51 ) becomes greater in error as f c increases since 
then Q c — > oo. 

An immediate solution to the problem is to substi- 
tute f c in place of Q c in eq . (fitf) , then obtaining 



A 



3/c 



mg - ^ Pi(r c ) 



(52) 



which gives the correct behavior of A e f f < nphom as 
t c — > for values of the clump filling factor. How- 
ever, now consider the behavior as t c — > oo, oc- 
curring when phom — > oo or k — > oo while hold- 
ing f c , a, and r c constant. Then Pi(r c ) — > 1 and 
Am g —> 3/ c /4r c as the clumps become opaque. For 
values of f c > 0.1 the resulting effective interaction 
coefficient A e // = 3/ c /4r c + npicm underestimates 
the actual value found from Monte Carlo simulations 
(see Figure |ll|(b) and discussion below). In addition, 
as f c — > 1 we should have A e jj — > nphom, but instead 
eq. (|5^) gives A e // — > 3/4r c + Kpi cm , which because 
of the dependence on r c will in general not reach the 
correct limit. 
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To extend the mega-grains approximation to all 
filling factors and retain the correct behavior as the 
clumps become optically thin or opaque, we propose 
to renormalize the clump radii by the factor (1 — / c ) 7 , 
where < 7 < 1 is a parameter that can fine tune 
the behavior as the clumps become opaque. In Ap- 
pendix |d| we discuss the motivation and a possible 
interpretation of this renormalization of the clump 
radii. Substituting r c (l — / c ) 7 for all instances of r c 
in eq.(p2h we get 



. _ 3^ P t [r c (l - / c ) 7 ] 



4rv 



(1 - UV 



(53) 



our new definition of the mega-grains interaction co- 
efficient. Now as t c — > or f c — > 1 we have that 
Pi [t c (1 — fc) 1 ] ~ | r c(l ~ /c) 7 and so the behavior of 
eq.(Ea) is 



A 



fc 



my 



\fic Picm) fc 



(54) 



since the (1 — / c ) 7 factor cancels out. This gives the 
desired result of A e // — > nphom as t c — > or / c — > 
1, and keeps A e /y < n phom for a clumpy medium, 
thereby extending the mega-grains approximation to 
the full range of < / c < 1. 

In the other extreme, when r c 3> (1 — / c ) -7 then 
-Pi [t c (1 - /c) 7 ] ~ 1 resulting in 



A, 



3/o 



4r c (l - / c ) 7 



(55) 



Comparing this with the previously discussed versions 
of the mega-grains approximation gives the following 
sequence of inequalities in the limit as the clumps 
become opaque: 



3/c 



4r c (l - f c ) 



3Q C 3/c . . 

> ~c > *r r ' (56) 



corresponding to equations (|53|), (ff6[), and ( (52|) re- 
spectively, with 7=1 and f c < 1. Varying the pa- 
rameter 7 allows for adjustment of the opaque clump 
limit behavior of eq.([53"|) over the complete range in 
the above sequence of inequalities: as 7 — > the left 
side of the inequality approaches the right side. We 
find that the optimal value of 7 in our extended mega- 
grains approximation depends on the arrangement of 
the source of photons with respect to the geometry 
of the clumpy medium: 7 ~ 0.5 gives a good ap- 
proximation for the situation of photons impacting a 
slab of clumps with a fixed angle with respect to the 



surface normal, whereas 7 ~ 0.75 works better for a 
central source within a sphere, and 7 = 1 works best 
for uniformly distributed sources. When / c < 1 the 
variation of the extended mega-grains approximation 
is small with respect to the parameter 7, and of course 
vanishes as the filling factor goes to zero. We shall use 
7 = 1 unless indicated otherwise. 

Figure ^ compares the three versions of the mega- 
grains approximation (MGA) to Monte Carlo sim- 
ulations. Figure |ll](a) shows the ratio T e ff/Thom 
versus Thom f° r the case of f c = 0.2, a — 100, 
and with r c = Q.05R, where R is the extent of the 
medium ( f c and a are dimensionless) . By definition, 
T e ff = RA eff and r hom = Rnp hom . The upper hor- 
izontal axis is the clump optical depth, r c , which is 
proportional to Thom via 

(a-1) 



(!) 



The 



(a - 1) f c + 1 



(57) 



derived by combining equations (|8|) and (|33J). The 
dashed line in Figure [ll](a) results from the use of 
equation (|4^) , the straightforward implementation of 
HP93, and it clearly exceeds unity for r c < 0.2 and 
Thom < 1, violating the requirement that r e // < Thom- 
The overshoot of unity becomes worse for larger fill- 
ing factors. The diamonds are results from Monte 
Carlo simulations of a central source in a sphere of 
radius R = 1 containing a two-phase clumpy medium 
characterized by the same parameters but with no 
scattering (ojd = 0). Each diamond plotted is the 
ratio T e ff I Thom, where r e // is computed by eq.(13). 
The dotted line results from the use of equation (52), 
giving T e ff < Thom as required, however the disagree- 
ment with Monte Carlo results becomes worse as Thom 
increases. The solid line results from the use of equa- 
tion ( |53| ) , our new definition of the mega-grains inter- 
action coefficient, also giving r e // < Thom and agree- 
ing well with the Monte Carlo results. 

Figure ^(b) compares r e ff resulting from the 
three version of the MGA at large T; lom (equivalent 
to large r c ) for the same clumpy medium. As before, 
the diamonds are results from the same Monte Carlo 
simulations. The dotted line is produced by using 
equation (p2h, giving 



Teff -> i?(3/ c /4r c + np lcm ) 



(58) 



as t c — > 00, which clearly disagrees with the Monte 
Carlo results. The dashed line is produced by using 
equation (E6h, giving 



T eff 



R(3Qc/4r c 



KPi, 



(59) 
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Fig. 11. — Comparison of three versions of the mega- 
grains approximations (MGA) with Monte Carlo re- 
sults, for the case of no scattering, in a clumpy 
medium having f c = 0.2, a = 100, and r c = 0.05. 
In both panels (a) and (b) the dashed line results 
from the original version of MGA using eq.(p6[), the 



dotted line from the modified version using cq.(52) 



and the solid line is our new extended version of MGA 
using eq.(|5^). The diamonds are results from Monte 
Carlo simulations. Panel (a) shows the ratio of the 
effective optical depth of the clumpy medium, r e //, 
to the equivalent homogeneous optical depth, T^ om , 
as a function of Thom and the clump optical radii r c 
(upper horizontal axis). Panel (b) shows r e // as a 
function of larger values of Thom and r c . 



as t c — > oo, which is in better agreement with the 
Monte Carlo results. The solid line is produced by 
our extended MGA, using eq. (|53|) , and as t c — > oo 



T eff 



R 



3/c 



4r c (l - fc 



+ KPi, 



(60) 



which provides the best approximation of the Monte 
Carlo results. There is a smooth transition in the 
slope of the T e jf curves as t c increases, and when the 
clumps become optically thick [r c > 4) their contribu- 
tion to the absorption becomes a fixed quantity so the 
only further change is due to the ICM density. Thus 
all versions of the MGA become linear as r c 3> 1 and 
Thom — ► oo, basically having the slope 



dr. 



eff 



dpi. 



1 



dT ho 



dpho 



(a-l)/ c + l' 



(61) 



which for the parameters used for Figure 0(b) is 
a value of about The intercepts of the asymp- 
totic lines with the vertical axis are given by R times 
the first term in each of the above eqs. 

(H©, and 

these values are approximately the average number 
of clumps encountered along a random line of sight 
(see Appendix |b|) . 

Now consider the behavior of the MGA as the 
clump radii vanish and other parameters are held 
fixed. This is shown in Figure 12 for the case of 
f c = 0.2, a = 100, and T ho m = 1, where r e // is 
plotted versus r c and r c . The dashed line produced 
by using eq. ( p6| ) exceeds the value of Thom — 1 (thin 
long dashed line) when r c < 0.03 (r c < 0.2), because 
as r c — > we have r c — > 0, and by eqs.(E8|) and (B^) 



T eff 



Rk l(p c 



l )Qc 



Thus the same problematic behavior of T e ff > Thom 
that was shown in Figure 11(a) for Thom ~ > occurs 
as r c — > and Th om is held constant. The solid line is 
produced by our extended MGA, using eq. (^3|) , and 
clearly yields the expected behavior r e // < Th om for 
the full range of clump radii, since then as r c — > we 
have the exact limit of 



T eff ~^ Rk l(Pc — Picm)fc + Pi, 



Tho 



[see eq.(p4j)1. The diamonds are results from Monte 
Carlo simulations, in good agreement with our ex- 
tended MGA. The dotted line is produced by using 
eq.(|5^), which has the correct limit as r c — > but 
tends to underestimate T e ff at larger clump radii. 
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Fig. 12. — Comparison of HP93 version [eq.(46), 
dashed line], modified version feq.(^), dotted line], 
and new extended version |eq.(|53|), solid line] of the 
MGA as a function of the clump radius r c . Results of 
Monte Carlo simulations are plotted as diamonds, and 
the clumpy medium is defined by parameters f c = 0.2, 
a = 100, and r hom = 1. 
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Fig. 13. — Comparison of r e // (vertical axes) pre- 
dicted by the new extended version [eq.(53|) with 
7 = 1, solid line], modified version [eg. (pi), dotted 
line], and original version [eq.([l6]), dashed line] of the 
MGA, with Monte Carlo results (diamonds), over the 
full range of clump filling factors (0 < f c < 1 on hor- 
izontal axis), for cases T% om = 1 and r^ om = 10, as 
indicated. The medium consists of clumps with radii 
r c = 0.05 and density ratio a — 100. 
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Figure U compares the extended MGA to the orig- 
inal version over the full range of clump filling factors, 
< f c < 1? with a — 100 and r c = 0.05 held constant, 
for cases r^ om = 1 and Tf lom — 10, as indicated. In 
each case the horizontal axis is f c and the vertical axis 
is T e ff, the effective optical depth. The diamonds are 
the T e ff resulting from Monte Carlo simulations of a 
central source in a sphere of unit radius containing 
a two-phase clumpy medium with clump filling fac- 
tors in the range 0.01 < f c < 0.8, also with a = 100 
and r c = 0.05 held constant. Clearly the extended 
MGA [using eq.(f~~ 



|) in eq.(|35[)1, shown as the solid 
lines, gives the best agreement with the Monte Carlo 
results, and satisfies the condition r e ff < Th om for the 
full range of clump filling factors. The original ver- 
sion of the MGA [using eq.@)], shown as the dashed 
lines, starts diverging to infinity for f c > 0.1 in the 
case of Thom — 1, and diverges to infinity for f c > 0.5 
in the case of Thom, — 10. The dotted line is produced 
using eq. (|52|) , which underestimates the value of T e f / 
and fails to reach the correct limit of T e f / — > Thom as 
f c — ► 1, especially when Thom is large. The fluctua- 
tions in the Monte Carlo results are due to variations 
in particular realizations of clumps in a finite volume, 
and these fluctuations grow larger as f c — ► 1. 

The behavior of r e // versus f c always exhibits a 
minimum at a single value of f c that has a compli- 
cated dependence on the other parameters character- 
izing the clumpy medium. To study this behavior, 
Figure |lj shows the prediction of r e / / by the extended 
MGA (solid lines) and the contribution from each 
phase on the clumpy medium: the dashed lines are 
ICM component, Ti cm = Rupicm, and the dotted lines 
are the component of the effective optical depth due 



tO Clumps, T mg = RA mg , SO that T e ff = T mg + T lcm . 

The horizontal axis is now log 10 f c and the cases of 
Thom — R^Phom = 1 and 10 are shown, with a = 100, 
and r c — 0.05 (same as Figure [l3]). Some of the 
same Monte Carlo results are shown (diamonds) to 
verify the behavior. Starting at f c ~ 0, Ti cm de- 
creases rapidly as f c increases since pi cm , given by 
eq. (H) , is inversely related to f c and the total mass is 
held constant. In fact the shape of the Ti cm = Rnpicm 
curve depends only on a and f c so only the magnitude 
changes with the value of Thom'- the shape of the vari- 
ation of Ti cm versus f c is exactly the same for the two 
cases of Thom m Figure [l4| The opacity of the clumps 
is also inversely related to their filling factor, as seen 
in eq. (|57|) , and so r c is near maximum when f c ~ 0, 
and then r. 
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Fig. 14. — Analysis of r e // versus f c predicted by 
our new extended version of the MGA for the same 
cases shown in Figure |l3|. The dashed lines are T mg = 
RA mg [eq. (153Q1) the dotted lines are Ticm — R^Picm? 
and the solid lines are the sum r e // = T mg + Ticm- 



mg 



R(4f c /3r c ). Thus T mg increases lin- 
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early with f c until the clumps become optically thin 
and they fill up the medium, causing r mg — > Thom- 
The linear increase of T mg does not compensate for 
the rapid decrease of r icm with f c , the net effect being 
that r e f f decreases with increasing f c until the clumps 
start to fill the medium. Therefore T e ff goes through 
a minimum as a function of f c when the number of 
lines of sight that pass through the very low density 
ICM greatly exceeds the number of paths intersecting 
optically thick clumps. For / C <1, the magnitude of 
T mg is independent of Thom whereas the magnitude of 
Ti cm is directly proportional to Thom, so the minimum 
of T e ff occurs at a larger value of f c for the case of 
Thom — 10 than for the case of Th om = 1- 

In general, the minimum of r e // occurs at nearly 
the same value of f c for which r mg — Ti cm , and analy- 
sis of the equations shows why this is true. Assuming 
that / C <1 then 



Or, 



my 



R 



df c " \3r c 
Assuming that a>l then 



dfc 



= -Rnp ic , 

Ticm 



' mg 



a — 1 



(62) 



(a-l)/ c + l 



fc + a 

Now since r e / / = T mg + Ti cm we have that 

^Teff ^~mg T C m 



dfc 



fc fc + — 



(63) 



(64) 



The local minimum occurs when the partial derivative 
is zero so that 



Ot, 



eff 



Ofc 



= 



1 + — T I r, 



mg 



, (65) 



which is approximately the relationship between the 
minimum of T e f / and the intersection of the T mg and 
Ticm curves in Figure |l4|. Consequently, as a function 
of fc (with constant a 1), the total luminosity 
absorbed by a clumpy medium attains a minimum 
when the luminosities absorbed in the clumps and in 
the ICM are equal, and this is demonstrated over a 



large region of parameter space in §7.2. 

The previous discussion of r e ff has been for the 
case of no scattering, or for scattering and absorp- 
tion together considered as interaction with dust. To 
model scattering in a clumpy medium we need the 



effective scattering albedo, w e //, given by eq. (p7|) , 
which depends on cu c , the clump albedo given by 
eq.(p4|). We find that the clump radius renormal- 
ization technique introduced above also needs to be 
applied to the clump albedo formula in order to pro- 
duce the correct approximation for scattering in the 
clumpy medium as f c — > 1. The new equation is then 



= UJ d V U esc [T C {1 - fc) 



U>d\ 



(66) 



so that lo c — > u>d as f c — > 1, since then the escape 
probability becomes unity. This leads to the correct 
behavior of u> e ff ~ > W <J f° r as fc ~ > 1- F° r the same 
reasons, r c in eq.(^) is replaced by r c (l — fc) 1 to give 
the correct limits as f c — > 1 for g c , the clump asym- 
metry parameter, and g e ff, the effective scattering 
asymmetry parameter. 

In summary, we have extended and improved the 
formulas for the effective optical depth and albedo 
of a clumpy medium, and introduced an approxima- 
tion for the effective scattering asymmetry parame- 
ter. The validity of the extended mega-grains formu- 
las, with scattering and for other source distributions, 
is further demonstrated by comparison with Monte 
Carlo simulations in and 

4.3. Comparison with Earlier Theory 

Here we compare the extended mega-grains ap- 
proximation (MCA) for the effect ive optical depth 
of just clumps with the theory of Natta fc Panagia 



1984)| (hereafter NP84), when a is large so the ICM 
can be ignored. Figure |l5| compares the approxima- 
tions and quantities involved as function of log 10 / c , 
for the case of r; lom = 2, a = 1000, and r c = 0.05. 
The upper panel shows r c (solid line), the clump op- 
tical radii, and two extreme values bounding F c , the 
average number of clumps encountered along a ran- 
dom line of sight: 



FS = R 



30c 
Arc 



(67) 



is the dashed line that eventually diverges to infinity, 
where Q c is the porosity given by eq.(44), and the 
dotted line is 

Ft 



(68) 



which stays finite, so that F[ < F c < F^ (see Ap- 
pendix |b|). F c is also called the covering factor of 
the clumps, and the random lines of sight are per- 
pendicular to the plane of a slab of thickness R, or 
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Fig. 15. — Comparison our extended MGA with the 
effective optical depth of a clumpy medium derived 
by Natta fc Panagia (1984)| , for the case of r\- lom = 2, 
a = 1000, and r c = 0.05. The top panel shows r c 
(solid line), the clump optical radii, or upper and 
lower bounds for F c , the covering factor of clumps 
(dashed and dotted lines). The solid line in bot- 
tom graph is r mg = RA mg produced by eq.( |53|). The 
dashed and dotted lines are produced by eq.(|69|), us- 
ing t c and bounds on F c from the top panel. The 
dash-dotted line is produced by eq(|70|). 



radial in the case of a spherical medium of radius R. 
In the lower panel the solid line is r mg — RA mg us- 
ing eq. (p3|) , the extended MGA. Assuming that the 
number of clumps encountered along a random line 
of sight is Poisson distributed with mean number F c , 
NP84 derived an equation for the effective optical 
depth: 

(69) 



(l - e - 4 ^ 3 ) 



where |t c is the expected optical depth of a random 
path through a spherical clump. The result of using 
F c = F^} in eq. (^9|) is shown by the dashed line in the 
lower panel of Figure |l^, and f agrees with r mg (solid 
line) for f c < 0.1, corresponding to F c < 1, but f 
exceeds r^ om when f c > 0.3 and F c > 5, diverging to 
infinity. If instead F c — F/ is used in eq.(|69]) then f 
remains finite as shown by the dotted line. Applying 
the clump radii renormalization gives 



- _ pf / l-cxp[-4r e (l-/ c )/3] 

1 - fc 



(70) 



shown as the dash-dotted line, which is almost iden- 
tical to r mg . So the approach used by NP84 gives 
almost the same result for r mg as the extended MGA 
when the average optical depth of a sphere, |r c , is 



used. The quantity 1 



-4t/3 



happens to be a good 



approximation of Pi(t), with a maximum relative dif- 
ference of 5% overestimation occurring at r w 1.5 
when Pj (1.5) » 0.8. 

4.4. Escape and Absorption Probabilities for 
Clumpy Media 

The effective optical depth, r e //, effective scatter- 
ing albedo, w e //, and effective asymmetry parame- 
ter, g e ff, given by the mega-grains approximation re- 
duces the basic radiative transfer properties of a two- 
phase clumpy medium to an effectively homogeneous 
medium. Therefore we propose use T e f f, u) e ff, and 
geff, given by equations ( j53| , |35| , |36| ), (|37| , |40| ), and 
(^2|, [i"3| ) respectively, in the escape and absorption 
probability formulae that were developed for homo- 
geneous media to estimate the escaping or absorbed 
fractions of photons in clumpy media. In particular, 
for isotropic emission uniformly distributed within a 
sphere we use P" sc given by eqs.(|l9|) and (pp|), for a 
central point source we use V^sc given by eqs.([l6|-|l8|), 
and for a uniformly illuminating external source we 
use V^ hs given by eq.(^). The full set of equations is 
summarized as a convenient list in ffi. 
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We find these analytic approximations give results 
in reasonable agreement with Monte Carlo simula- 
tions, as is demonstrated in As suggested by 
HP93, we also conjecture that given an escape proba- 
bility function for a homogeneous medium having any 
geometry and any distribution of sources, the escape 
probability for a two-phase clumpy medium of same 
geometry and source distribution can be reasonably 
approximated by substituting the effective optical pa- 
rameters given by the mega-grains equations into the 
homogeneous escape probability function. 

4.5. The Fractions Absorbed in Each Phase 
of a Clumpy Medium 

The energy of the radiation absorbed by the dust 
is converted into heat and then re-radiated as in- 
frared (IR) photons. The spectrum of the IR emis- 
sion from the dust is dependent on the dust tem- 
perature, which is most cases is determined by the 
equilibrium between the energy absorbed and emit- 
ted. In this paper, we shall not deal with very small 
dust grains, which can undergo temperature fluctu- 
ations, an effect that will cause a wider distribution 
in dust temperature and enhanced IR emission in the 
~ 3 — 20/im range (e.g. Dwek et al. 1997). For the 
case of a homogeneous medium with uniformly dis- 
tributed emitters, most of the dust will be at the 
same temperature (accept at the boundary of the 
medium where the temperature will be lower) . When 
the uniformly distributed emission occurs in an in- 
homogeneous medium there will be a distribution of 
dust temperatures. In a two-phase clumpy medium, 
the clumps will most likely attain a different radia- 
tive equilibrium temperature than the ICM, since the 
clumps will absorb/radiate a different amount of en- 
ergy. Given the total amount of energy absorbed by 
the dust we need to know what fraction is absorbed 
in the clumps and what fraction in the ICM, to de- 
termine the approximate equilibrium temperature of 
each phase of the clumpy medium. 

In the mega-grains approximation it is the over- 
density of the clumps that enters into the formulation, 
in this way separating the clumps and the ICM while 
keeping the ICM continuous. However, to determine 
the fraction of photons absorbed by the clumps we 
must consider the clumps as being completely sepa- 
rate from the ICM, therefore considering the full den- 
sity of dust in the clumps, and so the ICM is not 
viewed as uniform and continuous but as having holes 
occupied by clumps. This approach will give the cor- 



rect absorbed fractions as a — p c / Picm — * 1- Thus we 
define the full clump optical depth as: 

r c a = np c r c = anp icm r c , (71) 

where k is the total absorption plus scattering coeffi- 
cient of the dust grains per unit mass, and r c is the 
clump radius. We then replace t c with t° in all the 
mega-grains approximation equations, and the super- 
script "a" shall indicate that all of the clump optical 
depth is being used in the following quantities, in par- 
ticular, is the result of using r° in eq. (|34|) . 

First we deal with the case of a point source in 
the center of a spherical clumpy medium, and assume 
that the point source is not in any clump. Let R be 
the radius of the sphere and define optical depths 



-0 

mg 





icni 



_ u \a 
" ll mg 

= R Kpicn 



(72) 
(73) 



due to the mega-grains (subscript "mg") and ICM, 
respectively, where the superscript "0" indicates that 
this is for no scattering. To model the effects of scat- 
tering we employ the analytical approximation for the 
effective optical radius of a homogeneous sphere with 
scattering, as given by equations flT6| ) and (17), defin- 
ing new optical depths as 



'mg 
7~icm 



= r ° s (l-^)*(^c), (74) 
= ^ cm (l-co d )^^, (75) 



for the mega-grains and ICM, respectively, where the 
function x( T , 9) is given by eq.(|l7). Note that we use 
clump scattering properties u>% and g% for r mff , and 
the usual dust scattering albedo and asymmetry pa- 
rameter, Ud and gd, for the ICM optical depth. Of 
the total photons absorbed by the medium, the frac- 
tion, A\ cml that gets absorbed by the ICM is then 
estimated to be 



+ (l-/c)r i( 



(76) 



where the superscript "c" indicates a central point 
source. The factor 1 — f c is introduced to get an 
approximation that gives the correct behavior as a — > 
1 or f c — > 1, because the volume occupied by the 
clumps should not be included in the ICM. Now if 
Pe SC is the fraction of photons that escapes the sphere 
then A\ cm {\ — Vg SC ) is the fraction absorbed in the 



ICM and 
by clumps. 



— Pesc) is the fraction absorbed 
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In the case of an internal uniform distribution of 
emitters we consider the emission occurring in the 
clumps and ICM separately, as follows. For the frac- 
tion 1 — f c of the photons that are emitted in the ICM 
we predict that A c icm will be absorbed by the ICM. 
For the fraction f c that are emitted in the clumps we 
predict that a fraction 7- > " sc (t", uj c i) [given by eqs. (|Tgj) 
and (p0|)] of those photons will escape the clumps, and 
then a fraction A c icm of them will be absorbed by the 
ICM. The total fraction of emitted photons that will 
be absorbed by the ICM is then 

A u lcm ^ A\ cm [(1 - U) + fc V u esc {r a c ,uj d )] , (77) 

where the superscript "u" indicates that this is for 
the case of uniformly distributed emission of photons. 
This gives Af cm < A\ crn because part of the emission 
occurs inside clumps that are more dense than the 
ICM. Defining F abs = 1 - V^ sc {r e ff, u eff ) to be the 
total fraction of photons absorbed in the medium then 
A^ cm F abs is the fraction absorbed in the ICM and 
(1 — Af cm )F a bs is the fraction absorbed by clumps. 

To estimate absorbed fractions for the case of a uni- 
formly illuminating external source we just take the 
average of the values for central source and uniform 
source: 

Ac 

__icrn ( 7g ) 



Uniform Distribution of Emitters: 



si: 



Au 

icm 



The estimate is motivated by Monte Carlo simula- 
tions which indicate that Af cm < Af cm < A c icm . 
If Vabsi^ff^eff) is the total fraction absorbed in 
the sphere then Af cm V^ bs is the fraction absorbed 
in the ICM and (1 — Af cm )V* bs is the fraction ab- 
sorbed by clumps. The absorbed fractions predicted 
by the above equations are compared with the results 
of Monte Carlo simulations in §[| 

5. SUMMARY OF EQUATIONS 

5.1. Escape and Absorption Probabilities 

In the following formulae, use (r/j 0m , aid, fly) for all 
occurrences of (r, w, g) if the medium is homogeneous, 
or if the medium is clumpy use (r e / f,w e ff, g e f f ) given 
by the mega-grains equations below. 



Central Point Sou rce : 
PL( T > U >3) 

x{T,g) 



exp [-rs(T,w,g)] 
l-i(l- e - r/2 ) (1-9)- 



(79) 



Pc(r) 



l-w[l-P e (r)] 



P e (r) = -P(r) 



Pi(r) 



2r 2 



2r 2 



(80) 
(81) 
' 2t (82) 



Uniformly Illuminating External Source: 



V x abs (r,co) P(r) [l-< sc (r, W )] 

4r (1 - u) w , A . , 

= 3 l J esc( T ^) yP 6 ) 

Recall that P" sc and P^L are valid (most accurate) 
when g — g*(r) [see eq.(21)], whereas P e (r) and Pi(r) 
are exact formulae. 

5.2. Mega-Grains Approximation 

Input Parameters: 

k = dust mass extinction coefficient 

LOd = scattering albedo of dust grains 

gd = scattering asymmetry parameter 

Phom — average mass density of the dust 

f c = filling factor of the clumps 

7 = MGA tuning parameter ( ~ 1) 

a = Pc/ Picm 

r c = radius of each spherical clump 
Rs — radius of spherical medium. 



Effective Optical Depth, r e // : 



Ph> 



(ICM density) 



A cff 
T cff 



(a - 1) f c + 1 
(a — 1) Pi cm r c n (clump optical radius 

3/c Pi [r c (l - UV] 



(1 - fcV 



4r c 

Rs Kff 



(84) 



where the formula for A mg utilizes Pi(r) given by 
eq.(|8^) and < 7 < 1. Use 7 = 1 for uniform in- 
ternal/external sources, 7 = 0.75 for central source, 
and 7 = 0.5 for photons impacting a slab. 
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Effective Albedo, w e // : 



Ueff 



u c A mg + <jj d npi, 



A 



eff 



where the clump albedo is 

lo c = tu d V u esc [r c (l - f c )\ u d ] , 
and V^ sc is the OLEP formula stated in eq.(|80|) above. 
Effective Scattering Asymmetry Parameter, g e ff : 



9eff 



Qc ^-mg H~ 3d KPicm 



eff 



where the clump asymmetry parameter is 

1 + cxp(-B/A) 



= g d -C(l- 



l + exp([r c (l-/ c )T-S] I A) 
A = 1.5 + 4g d + 2uj dy /g^exp(-5g d ) 
B = 2 - g d (l - g d ) - 2u) d g d 



C 



/ 2g d ~- 2u d g d (l - g d ) 



5.3. Fractions Absorbed in Each Phase 

Let P a bs — 1 — Pesc be the generic absorption prob- 
ability. Then the fraction of photons absorbed in 
the ICM is Ai cm P a i, s , and the fraction absorbed by 
clumps is (1 - A icm ) P abs , where A icm is given by one 
of the following formulas, corresponding to the type 
of source. See below for definitions of r mg and Ti C7n . 



Central Point Source: 



4' 



(1 - U)tu 



+ (1 - fc)T« 



Uniform Distribution of Emitters: 

■^icm = ^icm [(1 " /c) + J 'c.Pesc( T c ') u '<*)] 

Uniformly Illuminating External Source: 



A x 



AU _i_ AC 

icm ' icm 



(85) 



(86) 



(87) 



Component Optical Depths: 



nig 
.0 



Kp c r c = 

Rs 



anpi, 



The superscript "a" indicates that the full clump op- 
tical depth is being used, i.e. A° lg is the result of 
using t" in eq. ([§4]) . The superscript "0" indicates 
that the quantities are for no scattering. To include 
the effects of scattering in spherical geometry, apply 
equations (79): 



>mg 



. (1 - u} a ) X ^9'Bt) 
' mg \ c ) 



r° m (l-^^™^). 



5.4. The Case of a < 1 

If one desires to model randomly distributed spher- 
ical cavities (e.g. SNRs instead of clumps) that have 
a lower density than the ICM (so that a < 1), the 
mega-grains equations can be applied by redefining 



f'o 

a' 



l-/ c 
1 

a 



(88) 



and using these new inverted values in the above equa- 
tions to obtain (r e ff,ui e ff,g e ff) for the medium with 
cavities. However, when computing the fractions ab- 
sorbed in each phase of the medium we do not ap- 
ply the inversion transform of eq.([S8|). The inverted 



mega-grains approximation is demonstrated in §3.2 



COMPARISON OF ANALYTIC APPROX- 
IMATIONS WITH MONTE CARLO SIM- 
ULATIONS 



In 54.2 



we developed the extended mega-grains 
approximation (MGA) and compared it with Monte 
Carlo simulations, which verified the effective opti- 
cal depth predicted by the approximations for the 
case of no scattering and when emission is from a 
central point source in a spherical clumpy medium. 
We shall use the acronym MCRT when referring to 
Monte Carlo simulations of radiative transfer, and the 
acronym MGEP when referring to the analytic ap- 
proximations consisting of the mega-grains (MG) ap- 
proximations combined with escape/absorption prob- 
ability (EP) formulae, which were summarized In 
this section we present more detailed comparisons of 
the MGEP model with MCRT simulations when scat- 
tering is also involved, and for all three types of source 
distributions discussed in In the following we des- 
ignate the cases of a central isotropic source by "C" , 
uniformly distributed internal sources by "U" , and 
uniformly illuminating external source by "X" . 
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We shall compare the MGEP approximations with 
MCRT results over a wide range of parameters that 
define a two-phase clumpy medium in spherical ge- 
ometry, except that the clump radii will be held fixed 
at r c = 0.05 relative to the radius of the sphere 
Rs = 1- A larger value would be less realistic and 
cause the number of clumps to be less than needed for 
good statistics in the MCRT model, whereas spherical 
clumps with radii smaller than r c = 0.05 become dif- 
ficult to represent accurately on a grid. All the Monte 
Carlo simulations were performed using a grid of 100 3 
or 127 3 voxels to represent the clumpy medium, and 
followed at least 10 6 photons for each result. The 
quantities in the comparisons will be shown as a func- 
tion of either equivalent homogeneous optical depth, 
Thom, clump filling factor, / c , or clump to ICM den- 
sity ratio, a, in the following three subsections. 

6.1. Dependence on Optical Depth 

Recall that the equivalent homogeneous optical 
depth of extinction, r^ om = nphomRs, can vary due 
to changing wavelength or dust mass (Rs is held con- 
stant), and so is the major parameter involved in 
modeling the transfer of a spectrum of radiation in 
a clumpy medium. Figure ^ showed the behavior of 
the effective optical depth t$ [eq.(13)] versus Thom 
resulting from Monte Carlo simulations with scatter- 
ing (ujd — 0.6 and gd = 0.6) for three cases of the 
filling factor, f c = 0.1, 0.2, and 0.3, (squares, trian- 
gles, and diamonds, respectively) with a = 100. The 
dashed, dotted, and solid lines are produced by the 
extended mega-grains theory combined with eq.(|79|) 
to compute Ts(T e ff,uj e ff,g e ff), and the curves agree 
well with Monte Carlo results. As discussed in § [T^ , 
the transition in the variation of t$ as a function r^ om 
goes from clump dominated (steeper slope at low val- 
ues of Thom) to ICM dominated (less slope at large 
Thom) as the clumps become opaque, and the occur- 
rence of this transition can be estimated by eq.(57). 
Unless otherwise indicated, a — 100 in this section. 



In Figure |16| we study the effect of changing the 
dust scattering parameters for the case of f c = 0.2, 
with the upper panel showing the effective optical 
depth with or without scattering (rs or r e ff respec- 
tively) for source type C, and lower panel showing 
the escaping fraction of source type U, all versus Thom- 
The case of no scattering is indicated with squares for 
MCRT results and the dash-dot lines arc the MGEP 
model. The diamonds and solid lines show MCRT 
and MGEP results, respectively, for u>d = 0.8 and 
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Fig. 16. — Comparison of MGEP theory (curves) 
to MCRT results (symbols) as function of Thom with 
f c = 0.2 and a — 100, for 3 different sets of dust 
scattering parameters as indicated in the plot legend. 
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g d — (isotropic scattering); the triangles and long- 
dashed lines show MCRT and MGEP results, respec- 
tively, for uj d = 0.8 and gd — 0.8 (almost forward 
scattering) . The values of T e f f given by MGA (dash- 
dotted line) for no scattering are used in eq. ([79]) , along 
with values of ^ e ff and g e ff (not shown), to com- 
pute Tg(r e ff,u)(,ff,g e ff), the solid and long-dashed 
lines in the upper panel. The solid line in the lower 
panel is computed by using r e // and u> e ff in eq.(|80|) 
for source type U. The MCRT results show that, as 
expected, forward scattering dust lowers the effective 
optical depth and increases the escape probability as 
compared to isotropic scattering dust, especially for a 
central source of photons, and the analytical approx- 
imations of MGEP correctly predict this effect. In 
the case of a uniformly distributed source the depen- 
dence of the escaping fraction on the dust scattering 
asymmetry parameter gd is weak, and so it is not a 
problem that g d does not enter into the MGEP equa- 
tions for source type U. The analytical approxima- 
tions of MGEP agree well with the numerical results 
of MCRT at low optical depths, whereas at high op- 
tical depths there are some differences, especially for 
source type C, which are due to specific realizations 
of the clumpy medium in MCRT simulations, but the 
agreement is still acceptable. 

Figure [l^ compares in the upper panel the ab- 
sorbed fraction of photons from each type of source, 
and in the lower panel the fraction of the total ab- 
sorbed photons that are absorbed in just the ICM, 
as function of Thorn, with f c — 0.1, and with scat- 
tering parameters u>d = 0.6 and gd = 0.6. The star, 
diamond, and square symbols represent MCRT nu- 
merical results for source types C, U, and X respec- 
tively. The solid, dot-dashed, and dashed lines repre- 
sent MGEP predictions for source types C, U, and X 
respectively. The absorbed fractions are computed as 

1 - T ? esc( T eff,^eff,geff) [C, eq. 

l-?£cfo//,we//) [U.eqjp] 

Kbs( T eff^eff) [X, eq.(|fj] 

for each indicated source type. The fractions ab- 
sorbed by the ICM, A icmi are computed by equations 
(76), (|77j), and ( |78| ) for source types C, U, and X, re- 
spectively. The MGEP model agrees with the MCRT 
results, showing the same behavior for each source 
type. Figure [l^ compares the same quantities for 
f c = 0.3. The MGEP model predicts a total absorbed 
fraction (upper panel) for source type X greater than 
the MCRT results when Thom > 10, probably due to 
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Fig. 17. — Comparison of MGEP absorption prob- 
abilities to MCRT absorbed fractions as function of 
Thom, with f c — 0.1, a — 100, and r c = 0.05. The top 
panel shows the total absorbed fractions whereas the 
bottom panel shows the fraction of the total which 
is absorbed by the ICM. Scattering parameters are 
LUd = 0.6 and g& = 0.6. The star, diamond, and square 
symbols represent MCRT results for source types C, 
U, and X, respectively. The solid, dot-dashed, and 
dashed lines represent MGEP results for source types 
C, U, and X. 
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Fig. 18.— Same as Figure but for f c = 0.3. 



limitations of the assumptions used in deriving eq. (|29|) 
for V*bs- For this case of f c — 0.3 (and actually for 
f c > 0.1) the MGEP model predicts values for A icm 
(lower panel) that are lower than the MCRT results. 
When f c < 0.1 the absorbed ICM fractions predicted 
by MGEP are slightly higher than the MCRT results. 
Overall, the agreement is acceptable and the MGEP 
model exhibits the same behavior as the MCRT sim- 
ulations. 

6.2. Dependence on Filling Factor 

We next compare MGEP and MCRT results as a 
function of the clump filling factor f c , holding other 
parameters fixed. Figures |l9| and compare MCRT 
and MGEP predictions of the fraction of emitted pho- 
tons that are absorbed in the medium, P a b s , in the up- 
per panels, and the fraction of the total absorbed pho- 
tons that are absorbed in the ICM, Ai cm , in the lower 
panels, as function of f c , with scattering parameters 
uJd = 0.6 and gd = 0.6. Figure |l^ is for the case of 
Thom — 2 and Figure ^ is for Thom = 10 (note the dif- 
ferent ranges for P a ^ s axis). The star, diamond, and 
square symbols represent MCRT numerical results for 
source types C, U, and X, respectively. The solid, dot- 
dashed, and dashed lines represent MGEP analytical 
results for source types C, U, and X, respectively. The 
analytical MGEP theory agrees well with the MCRT 
numerical results, with the exception of the absorbed 
fraction of a uniformly illuminating external source 
with scattering. For that case of source type X, when 
Thom is large and f c — * 1 the MGEP theory overes- 
timates the absorbed fraction with a difference that 
increases with f c . This is probably because the as- 
sumption of uniformly distributed photons after first 
scattering (see §3.4) becomes invalid as Thom 00 
since then the externally impacting photons do not 
penetrate the sphere. 

The dotted line and the triangles in the lower pan- 
els represent MGEP and MCRT results, respectively, 
for the case of no scattering {ojd = 0). The effect 
of turning on scattering is to decrease the fraction of 
photons absorbed in the ICM, since photons scattered 
by the ICM then have more of a chance of getting ab- 
sorbed in the denser clumps. The MGEP model tends 
to overestimate A icm at low clump filling factors and 
underestimate Ai cm at high filling factors, with best 
agreement at about f c — 0.1. However, the MGEP 
theory gives the same variation of A; cm with f c as 
the MCRT simulations, and overall the agreement is 
acceptable. 
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Fig. 19. — Comparison of MGEP absorption proba- 
bilities to MCRT absorbed fractions as function of f c , 
with Thom = 2, a = 100, and r c = 0.05. Scattering 
parameters are LUd — 0.6 and gd — 0.6. The star, di- 
amond, and square symbols represent MCRT results 
for source types C, U, and X, respectively. The solid, 
dot-dashed, and dashed lines represent MGEP results 
for source types C, U, and X. The dotted line and the 
triangles represent MGEP and MCRT results, respec- 
tively, for the case of no scattering. 
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Fig. 20. — Same as Figure [La but for Thom = 10. 
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6.3. Dependence on Density Ratio 

Finally, we compare MGEP and MCRT as a func- 
tion of the clump to ICM density ratio, a = p c / Picmi 
holding other parameters fixed. Figure ^l] shows the 
effective optical depth (r e // or rs) for source type 
C (top row), the escaping fraction of photons (Pe SC ) 
from source type U (second row), the fractions ab- 
sorbed by the ICM for source type C (Aj cm in third 
row), and for source type U {A" cm in fourth row), all 
as a function of a, for the cases of f c = 0.1, 0.5, and 
0.9 as indicated, with r^ om = 10. The square and 
diamond symbols indicate MCRT results for no scat- 
tering (ujd — 0) and with scattering ( for u>d = 0.6 
and gd = 0.6), respectively. The solid and dashed 
lines indicate MGEP predictions with and without 
scattering, respectively. The value a = 1 corresponds 
to a homogeneous medium, and then r e ff = Th om 
when LOd — 0. Values of a < 1 show the case of spher- 
ical cavities in a denser ICM. For the case of denser 
clumps (a > 1), most of the variation in r e //, rs, and 
■p" sc occurs for 1 < a < 100 because in this range 
the clumps go from being optically thin to thick, and 
for a > 100 the clumps are essentially opaque. Turn- 
ing on scattering decreases the effective optical depth 
and allows more photons to escape, but this effect is 
diminished as a — > oo, and is correctly predicted by 
the MGEP model. The reason for the diminished ef- 
fects of scattering is that ui e ff and <? e // (not shown) 
decrease by more than a factor of two as a — ► oo. 

Applying the MGEP model to the case of spheri- 
cal cavities (a < 1) is accomplished by the inversion 
transform of eq.(^8(), which then regards the cavities 
as the ICM and the rest of the medium as clumps. 
This is not the intended application of the mega- 
grains approximation, since there are then actually 
no isolated clumps, but there is some indication that 
the MGEP model could be tailored to give a reason- 
able approximation of the escaping radiation. When 
f c = 0.9 we see that the MGEP model does not give 
correct predictions for a < 1, the case of spherical 
cavities. This is understandable, because when the 
cavities occupy 90% of the volume, the high density 
regions are just shells between the cavities, which are 
not well approximated by clumps having the same 
radius as the cavities (recall r c = 0.05). Choosing 
a smaller value for r c would improve the MGEP ap- 
proximation of the MCRT results. The MGEP model 
provides an acceptable approximation of the MCRT 
results for a > 1. When f c = 0.5 and as a — > oo there 
is an increasing difference between MGEP and MCRT 



results in the case of a central source, but this can be 
attributed to the specific realization of the clump lo- 
cations relative to the source: there is a 50% chance 
that the point source will be inside a clump (as in this 
case) , causing an increase of the optical depth in the 
MCRT simulation. Other simulations for the same 
filling factor give MCRT results that are less than the 
MGEP predictions, so we expect that on average the 
MGEP predictions will be close to MCRT results. 

The fractions absorbed by the ICM are shown in 
the third and bottom rows of Figure ^ for source 
types C and U, respectively. The intersection of 
the horizontal and vertical dotted lines indicate the 
nominal value of Af cm expected as the medium be- 
comes homogeneous, since then the fraction photons 
absorbed by the ICM should be 1— f c of the total ab- 
sorbed photons. For the case of when the clumps are 
denser than the ICM (a > 1) the MGEP model agrees 
well with the MCRT results, when the source is uni- 
formly distributed. When the source is at a central 
point and f c = 0.5 the MGEP predictions are less 
than the MCRT results, and this can be attributed 
to the fact that the point source happens to be in 
a clump, as discussed above. Note that the MGEP 
model correctly predicts that turning on scattering 
causes less photons to be absorbed by the ICM and 
more to be absorbed in clumps. When a < 1 the 
MGEP model predicts values lower than MCRT for 
f c > 0.3, since the MGA was not actually formulated 
for this application, but the MGEP model does ex- 
hibit the correct behavior. 

6.4. Body Centered Cubic Lattice of Clumps 

The mega-grains approximation can be applied to 
the case of cubic clumps randomly located on a body 
centered cubic (BCC) lattic e with partial success. 
The BCC lattice was used by |Witt fc Gordon (1996) 



in their Monte Carlo simulations. For source type 
U, the MGEP model does predict the same escaping 
fraction computed by MCRT when r c — 1/N = 0.05, 
where N = 20 is the number of grid element along 
each axis of the lattice. However, this value of r c 
in the MGEP model gives erroneous predictions for 
the effective optical depth seen by a central source. 
The problem is related to the fact that the appar- 
ent projected area of a cubic clump depends on the 
viewing angle. To match the MGEP model with 
MCRT results for source type C, we need to use 
r c ~ \/2/N — 0.071 in the mega-grains equations, 
but then the predictions for source type U are wrong. 
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Fig. 21. — Comparison of MGEP and MCRT results as f unc tion of a = p c /picm for three values of f c = 0.1, 0.5, 
and 0.9, as indicated, with Thom — 10 and r c = 0.05. See §6.3 for more information. 
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7. MODELING THE ABSORPTION OF 
STELLAR RADIATION BY DUST AND 
THE INFRARED EMISSION 

In this section we model the transfer of a spectrum 
of radiation in a two-phase clumpy medium, from 
emission by starlike sources to the absorption and 
scattering by dust and the resulting infrared emission 
from dust heated by the absorbed radiation. Both de- 
tailed Monte Carlo simulations and the analytical ap- 
proximations developed in previous sections are used 
to model the transfer of radiation, and the emerg- 
ing spectral energy distribution (SED) from the two 
methods of modeling are compared. We shall again 
use the acronyms MCRT when referring to the Monte 
Carlo radiative transfer model, and MGEP when re- 
ferring to the mega-grains approximations combined 
with escape/absorption probability formulae. 

The dust is assumed to be composed of 40% graphite 
and 60% silicates by mass, and to have the grain size 
distribution 

C(o) oc a" 3 ' 5 , (89) 
normalized over the following range of grain sizes: 

0.001/xm < a < 0.25^m . (90) 



Using the optical constants from Draine (1985), we 
applied Mie theory to calculate the absorption and 
scattering efficiencies of dust grains, Qabsifl, A) and 
Qscat(a, A) respectively, as a function of grain size, a, 
and photon wavelength A. The scattering asymme- 
try parameters, g(a, A) = {cosO sca t(a 1 A) ), were also 
calculated by averaging with respect to the distribu- 
tion of scattering angles. Then the cross-sections and 
asymmetry parameters were averaged with respect to 
the grain size distribution to get the mass absorption 
and scattering coefficients and average asymmetry pa- 
rameters used in the models: 



{Tra 2 Q abs (a, A)} 

KabsW = -, r — 

( m 9>a 

(Tra 2 Q scat (a,X)) a 



(91) 



K SC at{X) 



9dW 



_ (g(a, \)ira 2 Q scat (a, A)) q 



(na 2 Q scat (a. X)) a 
where the averaging operator is defined as 



(■>«. 



(•) ((a) da 



(92) 



and the average grain mass is 

4 



ira 3 p 



(93) 



where p is the mass density of a graphite or silicate 
grain. All the following simulations have the same 
total dust mass of 1.01 Mq contained in spherical re- 
gion lpc in radius. This corresponds to a dust mass 
density of ptom — 1-63 x 10~ 23 gm/cm 3 in the homo- 
geneous case, equivalent to a gas density of 1000 cm~ 3 
with a dust to gas mass ratio of 0.007, giving a ho- 
mogeneous optical depth of ry = 1.67 in the V-band. 

Before discussing the results of the MCRT and 
MGEP simulations, we illustrate how the degree of 
dumpiness in the distribution of dust is as important 
as the total dust mass and scattering albedo in af- 
fecting the transfer of radiation through the medium. 
Figure ^2| compares as a function of wavelength, (a) 
the effective optical depth, r e //(A), (b) the effective 
albedo, u; e //(A), and (c) the effective asymmetry pa- 
rameter, g e ff(X), of a spherical region of dust with 
different degrees of dumpiness, computed using the 
mega-grains equations listed in §[| with Rs = lpc, 
= 1.63 x 10^ 23 gm/cm 3 , and 



Ph 



«(A) = KabsW + K scat (X) 

^d(A) = K scat (X) / k(X) 



(94) 



The solid curves in Figure ^2] show the homogeneous 
case, giving ry = 1.67. The dashed curves are for the 
case f c = 0.1 with a — p c /picm = 100, and clump 
radii r c = 0.05pc, giving ry = 1.11. The dotted 
curves are for the extreme case of f c — 0.01 with 
a = 10 4 , and r c = O.Olpc, giving ry = 0.66. It is 
evident that the effective radiative transfer proper- 
ties of the dusty medium can be radically affected by 
the degree of dumpiness. The effect is greatest at 
shorter wavelengths where the dust absorption and 
scattering coefficient is sufficiently large to make the 
clumps opaque, causing r e //(A) to become almost flat 
and featureless. The small features remaining are due 
to the ICM, but the opaque clumps dominate the ef- 
fective optical depth. An explanation for this "gray" 
behavior when r c (A) ^> 1 is obtained by using eq. (|6l|) : 

dT eff 



OX 



dl ~eff 






dX 








dX 



(95) 



dX 



Thus, as af c — > oo the variation of r e //(A) for 
the clumpy medium becomes a small fraction of the 
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Fig. 22. — Comparison of radiative transfer proper- 
ties of dust distributed homogeneously or in clumps. 
Panels show effective values of (a) optical depth, (b) 
albedo, and (c) scattering asymmetry parameter, ver- 
sus wavelength, computed using the mega-grains ap- 
proximation. For comparison, the solid curves show 
the dust optical parameters for the equivalent homo- 
geneous medium. The dashed curves show effective 
optical parameters for a two-phase clumpy medium 
with filling factor f c — 0.1, clump radii r c — 0.05pc 
(843 clumps), and clump to ICM density ratio of 
a = Pc/picm — 10 2 . The dotted curves show the 
case of f c — 0.01, r c = O.Olpc (9918 clumps), and 
a = 10 4 . See §7 for more information. 



variation of T/ lom (A) in the equivalent homogeneous 
medium. For longer wavelengths (infrared) there is 
no difference between clumpy and homogeneous me- 
dia because the clumps are optically thin. 

To test the MGEP model using MCRT, we sim- 
ulated a two-phase clumpy medium with f c — 0.1, 
a = 100, and r c = 0.05pc, so that the effective opti- 
cal depths, albedos, and asymmetry parameters used 
in all the MGEP calculations are given by the dashed 
curves in Figure |2^. The MCRT computations always 
use the optical parameters given by the homogeneous 
case (solid curves) in Figure [2^, simulating the details 
of the clumpy medium on a 3D grid of 127 3 voxels. 
The radiation source in all cases is a black-body spec- 
trum with T s = 15000K and L s = 33000L Q , so the 
spectrum of the emitted flux is 

Sx=-^B X (T S ). (96) 

The transfer of radiation is computed by MCRT 
at 40 wavelengths from A m ; n = 0.1 (im to A max = 
14/zm, following 10 7 photons at each wavelength. The 
standard three types of source geometries are stud- 
ied: uniformly distributed internal emission (U), uni- 
formly illuminating external source (X) , and a central 
isotropic point source (C). 

7.1. Absorbed Luminosities and the 
Distribution of Dust Temperatures 

Since we model only non-ionizing photons experi- 
encing coherent scattering, for a given type of source 
the radiative transfer of a unit emitted flux can be 
simulated at each wavelength separately, giving an 
escape and absorption response function for the par- 
ticular choice of source, geometry, and dumpiness. 
Then wc multiply the chosen source spectrum with 
the escape/absorption response function to get the 
actual escaping SED and the luminosity absorbed by 
each component of the dust. In the case of Monte 
Carlo simulations (MCRT) we obtain the 3D spatial 
distribution of the absorbed luminosity, whereas us- 
ing the analytical approximations (MGEP) gives the 
luminosity absorbed by all the clumps and luminosity 
absorbed by the ICM. The equilibrium dust temper- 
atures, Td for graphite or silicate, are computed by 
equating the absorbed and emitted luminosities: 

/•A max /.oc 

/ S x P abs (X)d\ = 4Trm d Kx B x (T d )dX, 

(97) 
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where P a bs(A) generically represents the fraction of 
flux at wavelength A absorbed by graphite or silicate 
in a particular 3D voxel of mass in the case of 
MCRT, or it is the fraction absorbed by graphite or 
silicate in either the ICM or the clumps of mass in 
the case of MGEP. The energy conservation eq. ( |97j ) is 
solved iteratively for Td using Brent's method of find- 



ing the zeros of a function (Press et al. 1992, pp.242 
254])], giving a 3D distribution of dust temperatures 
for MCRT, or the average temperatures of the dust 
in clumps and the ICM when using the MGEP model. 

Figure |23| shows the distribution of equilibrium 
dust temperatures computed by the MCRT model for 
clumps (top row) and in the ICM (bottom row), for 
source types U, X, and C in each column from left 
to right. The solid lines are the probability densities 
of graphite dust temperatures and the dotted lines 
represent silicate dust. For the cases of uniformly 
distributed internal and external sources the proba- 
bility densities are Gaussian on the high temperature 
side and of exponential form on the low temperature 
side (note slightly different temperature scales). In 
contrast, the central source creates a power law dis- 
tribution of temperatures (third column is a log-log 
plot) since dust next to the source is heated much 
more than dust farther from the center. In fact, the 
straight lines in Figure ^3] are power laws: the solid 
line is TJ 8 and appears to be parallel to the graphite 
temperature distribution whereas T~[ 9 (dashed line) 
is parallel to the temperature distribution of silicates. 

Table | compares the luminosities absorbed by 
each component (graphite or silicates) and phase 
(ICM or clumps) of the medium, as computed by the 
MCRT and MGEP models, again for the cases of uni- 
form internal, uniform external, and central source. 
The MGEP results were obtained by computing the 
escaping and absorbed SED using the mega-grains 
equations and the escape/absorption probability for- 
mulae, combined with formulae for the fra ctio ns ab- 
sorbed by the ICM and clumps, given in §|4.5|. Also 



compared in Table are the resulting dust tempera- 
tures: averages of the distributions shown in Figure |23| 
for the case of MCRT simulations, and in the case of 
the MGEP model when the source is not a central 
point source, a single temperature for each compo- 
nent and phase computed using eq.(|9^). 

To compute the average temperature of the dust 
in the case of a central source we use a theoretical 
power law distribution of dust temperatures as de- 
scribed in Appendix [e], along with two assumptions: 



that the absorbed luminosity decays with radial dis- 
tance from the source like an inverse power law with 
index 77 = 2.5 (77 = 2 is normal for optically thin 
case), and we use the maximum dust temperatures 
found by MCRT in the voxel containing the source. 
Normally the maximum temperature would be the 
dust sublimation temperature, however, the MCRT 
simulation is not set up to fully resolve the volume 
next to the source so the sublimation temperature is 
not reached in the simulation. Since we are compar- 
ing with MCRT results we chose to use those maxi- 
mum temperatures. We then solve for minimum tem- 
perature of the power law distribution that matches 
the emitted and absorbed luminosities, and this also 
yields the average dust temperature. 

Table shows that the MCRT and MGEP re- 
sults are in close agreement. The uniformly illumi- 
nating external source experiences more absorption 
(67%) than the uniformly distributed internal emis- 
sion (49%) and therefore heats the dust to slightly 
higher temperatures. Approximately the same frac- 
tions of luminosity are absorbed in the central and 
external source cases, but the average dust tempera- 
ture is lower in the central source case since most of 
the dust is far from the source. In all cases the clumps 
absorb about seven times more energy than the ICM, 
and this is in part due to the fact that the clumps 
contain more than ten times the mass of the ICM: 



M r 



PcfcV 



(1 - fc)V 



1-/0 



(98) 



and this ratio is greater than ten for f c = 0.1 and 
a = 100. However, the larger mass can easily radiate 
the absorbed energy in the IR so that the temperature 
of dust in clumps is generally lower than in the ICM. 
In §7.3 we study in more detail the reasons for the 



lower temperatures of dust in clumps. 
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Table 1 

Comparison of Dust Energetics and Temperatures 









L abs (1O 3 L0) 


Tdust 


(K) 


source 


dust 


phase 


MCRT 


MGEP 


MCRT 


MGEP 


U 


graphite 


ICM 


1.66 


1.56 


45.8 


45.7 






Clumps 


11.55 


11.54 


41.9 


42.2 




silicates 


ICM 


0.44 


0.42 


37.0 


37.0 






Clumps 


2.75 


2.74 


33.4 


33.6 




total 


Labs/ L s = 


49% 


49% 






X 


graphite 


ICM 


2.33 


2.19 


48.7 


48.5 






Clumps 


15.70 


15.86 


44.2 


44.6 




silicates 


ICM 


0.62 


0.58 


39.0 


39.0 






Clumps 


3.66 


3.69 


34.9 


35.3 




total 


L a bs/L s = 


67% 


67% 






c 


graphite 


ICM 


2.40 


2.21 


41.3 


42.8 






Clumps 


15.16 


15.48 


37.9 


39.9 




silicates 


ICM 


0.64 


0.59 


32.4 


34.8 






Clumps 


3.57 


3.63 


29.3 


31.9 




total 


LabslL s = 


65% 


66% 
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Fig. 23. — Probability distribution of equilibrium dust temperatures from MCRT simulations for each type of 
source (indicated at top), in a clumpy medium with f c = 0.1, a = 100, and r c = 0.05. The top row represents the 
temperature distribution of dust in clumps and the bottom row that for the ICM. The solid lines are for graphite 
and the dotted lines are for silicates. See £7.1 for more information. 
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Fig. 24. — Comparison of dust emission spectra from MCRT and MGEP models for each type of source: U, X, C, 
from left to right (indicated at top). The clumpy medium has f c — 0.1, a = 100, and r c — 0.05. The top row is 
for graphite, bottom row is for silicates. The triangles are for the MCRT model of dust in ICM and the diamonds 
are for dust in clumps. The solid lines are for the MGEP model of dust in ICM and the dotted lines are for dust 



in clumps. See §7.2 for more information. 
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7.2. The Emerging UV-FIR SED 

Computation of the dust temperatures also gives 
the IR emission from dust in each voxel of the 3D 
simulation by the MCRT model, or from each phase 
of the medium in the MGEP model: 

Fx = inm d K(X)B x (T d ) , (99) 

where T d is the temperature of dust in a voxel of mass 
ma in the MCRT model, or where T d is the temper- 
ature of a dust in clumps/ICM of mass rrid in the 
MGEP model. Equation @ is used by the MGEP 
model only for uniformly distributed internal or ex- 
ternal sources. In the case of a central source with a 
power law distribution, p(Td), of dust temperatures, 
the MGEP model essentially computes 

Fx = A-Km d K{\) I ^ B x {T)p{T)dT, (100) 

JT min 

an approximation that is described in more detail in 
Appendix [e]. Integrating the IR emission over the 3D 
volume in the case of the MCRT approach, or just 
adding the IR emission from the clumps and ICM in 
the case of the MGEP approach, gives the total IR 
emission spectrum. Figure |^ shows the IR emission 
spectra from graphite (top row) and silicates (lower 
row), with the diamonds and triangles representing 
emission from clumps and the ICM, respectively, as 
computed by MCRT, and the solid and dotted lines 
representing emission from the ICM and clumps as 
computed by MGEP. The source types U, X, and 
C, are presented in columns from left to right. The 
MGEP model is in close agreement with the MCRT 
results. The emission from the clumps is in general 
greater than that from the ICM because the absorbed 
luminosity is larger. An exception is the case of a 
central source where the heating of dust adjacent to 
the source to much higher temperatures causes more 
emission from the ICM at short IR wavelengths. 

The emerging SED is the sum of the escaping ra- 
diation and the IR emission from heated dust: 

XE X = XSxPescW + AF A , (101) 

where P esc (X) generically designates an escape prob- 
ability that is computed either numerically (MCRT) 
or analytically (MGEP) and depends on the model 
and source type. Figure ^5] compares the SED com- 
puted by the MCRT and MGEP models for each type 
of source. The diamond symbols are the MCRT re- 
sults, the solid lines are the MGEP theory, and the 



XE X (ergs/s) 




0.1 1.0 10.0 100.0 

X (yU,m) 



Fig. 25. — Comparison of the emerging UV-FIR SED 
resulting from the MGEP model (solid lines) and the 
MCRT model (diamonds). The dashed line is the 
SED of the source types U, X, and C, from top to 
bottom rows. 
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dashed line is the original source spectrum XS\ . The 
source types U, X, and C, are presented in panels from 
top to bottom. The SED at short wavelengths is the 
source radiation that escapes from the clumpy dusty 
environment, whereas the SED at long wavelengths is 
emission from dust heated by the absorbed radiation. 
The agreement between the MGEP and MCRT mod- 
els is excellent. The SEDs for the uniform internal 
and external sources have very similar appearances, 
except that in the case of an external source more of 
the source energy is absorbed (68% versus 49%). The 
central and external source cases have the same frac- 
tion of energy absorbed by the dust, so that the short 
wavelength portion of the SEDs are the same. How- 
ever, in the case of a central source the emission from 
hot dust near the source increases the SED around 
10/xm. In all cases, the relative paucity of emission in 
the ~ 3 — 20/um region is caused by the omission of 
stochastic heating in the calculations. 

7.3. Exploring Parameter Space with the 
MGEP Model 

Since the MGEP equations are computationally 
fast, we can very easily model the escape and ab- 
sorption of radiation over a wide range of parame- 
ters characterizing the clumpy medium or the radi- 
ation source. Figures and |27] show the results of 
MGEP model calculations over the range of clump 
filling factors, 0.001 < f c < 1, and density ratios, 
1 < a < 10 4 , with clump radii r c — 0.05. The 
dust composition and total mass is the same as that 
used in the above MCRT to MGEP comparisons. The 
sources of radiation are isotropic and uniformly dis- 
tributed within the sphere, again having a black-body 
spectrum with temperature T s = 15, 000K and total 
luminosity L s — 33, 000-L©. In all the panels the dia- 
mond and cross marks the point (a, f c ) — (100, 0.1) 
for which we presented detailed comparisons with 
MCRT, and the horizontal dotted lines indicate the 
limiting value of f c (=0.0015) below which there are 
less than 10 clumps in the medium, therefore, only 
filling factors above the dotted line are considered to 
be statistically reliable. 

Figure ^ shows (upper panel) the contours of the 
fraction of the source luminosity that is absorbed by 
the clumpy medium (L a t> s /L s ), and (lower panel) the 
contours of the logarithm of the ratio of luminosity 
absorbed by the ICM to that absorbed by clumps 
[\og(Li cm / L c )]. The dashed line (in all panels), ex- 
tending from (aj c ) = (1,0.5) to (1000,0.001), is 
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Fig. 26. — Results of the MGEP model over a wide 
range of clumpy medium parameters (a, f c ) for a uni- 
formly distributed black-body source having T s = 
15, 000K and L s = 33, 000L©. (a) Upper panel shows 
absorbed luminosity fraction; (b) and the lower panel 
shows the ratio of ICM/clump luminosities. The dia- 
mond and cross marks the parameters for which de- 
tailed comparisons with MCRT were presented. In 
both panels the dashed line is the locus of parameters 
where M icm = M c , and the thick solid line is the lo- 
cus of pa ram eters where Li cm — L c . More details are 
given in §7.3 of the text. 
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the curve 



/m(«) 



1 



(102) 



giving the locus parameters for which the ICM and 
clumps have equal mass, that is, M icm = M c . The 
equation for /m(o0 is obtained from eq. (|98|) , which 
we restate here as 



M n = 



l-/c 



M,, 



(103) 



Note that if f c > f M {a) then M c > M lcm . The thick 
solid line in all panels represents the locus of parame- 
ters /l(«) for which the ICM and clumps absorb the 
same amount of energy, that is, Li cm = L c . The dot- 
ted contours in the lower panel of Figure indicate 
when the ICM absorbs less energy than the clumps. 
The solid contours in the lower panel of Figure [2^ in- 
dicate that if f c < fh (a) then L c < Li cm . Combining 
this with the fact that M c > Mi cm when f c > /m(o;) 
we have that for the parameter region between the 
thick solid and dashed curves 



fin (a) < fc < h{a) 



M c > M icm 



(104) 



which guarantees that the temperature of dust in 
clumps will be less than that in the ICM since, com- 
pared to the ICM, a smaller amount of energy is ab- 
sorbed by a larger number of dust particles, which 
can radiate this energy at a lower temperature. 

Figure (upper panel) also shows that for a given 
density contrast ao, the minimum absorbed luminos- 
ity {L a bs) corresponds to the filling factor f c = /z,(ao) 
for which L c = Li cm . When a ^> 1, the condi- 
tion L c = Li cm is achieved when there is approxi- 
mate equality between the effective optical depth of 
the clumps and that of the ICM, that is, T mg (X p ) » 
Ticm{\), where X p is the wavelength where most of 
the source energy is absorbed. By equation (|65|), the 
value of f c for which r mg « Ti cm is nearly the same 
value for which r e / / attains a minimum as a function 
of f c (see also Figures [l4|, [ll] and [2(]) . Consequently, 
as a — > oo, the steepest descent of L aos is along the 
curve /z,(a). A detailed derivation of this result and 
the following equations would require studying T mg 
and T icm as given by eqs.(Q) and (f75|), respectively, 
but here we assume that the effects of scattering are 
of second order importance. 

Let us study the curve /l(q;) in more detail. As 
a — > 1 the condition L c = Li cm is achieved by having 
M c = Micm, and so ~ fni(a) when a < 5, as 



shown by the merging of the dashed and thick curves 
in in Figure However, as a — ► oo we see that 
fh(ct) ^> /m(o)- This behavior is due to the change 
in the effective optical depth of the clumps, r mg , since 
the clumps become optically thick as a — + oo, so we 
derive an equation for /l(«) as follows. Starting with 
T m g(\) = T icm (A )9 ) and applying the mega-grains 
equations from we obtain 



3/c 



K(\p)Pha 



4r c (l - foV {oc - 1) fc + 1 



(105) 



assuming optically thick clumps. Since we propose 
that eq.(105) occurs when L c — Li cm we now use 
the symbol /l in place of f c for the solution. Then 
as a — > oo we know that /l <C 1 but afi 3> 1, so 



eq.(lOS) becomes 



3/l 

4r c 

Solving for /x, gives 



fUa) 



k(X p )p, 



on i 



Oifl 



4r c K,(X p )p h 

o 



3o 



(106) 



(107) 



when a is large. Using K(X p )phom — 1-7, which occurs 
when X p ~ 0.55/im, gives fh{ot) ~ 5 Q ^ which pro- 
duces a perfect fit to the zero contour of log(Lj cm / L c ) 
in Figure ^6| (lower panel) when a > 30 (not shown). 
The transition from the behavior /i(a) ~ fAi(ct) to 
the behavior /i(a) a occurs naturally around 
values of a for which eq.( |l07[ ) is greater than /m(o;), 
indicating that the clumps are becoming optically 
thick. We find that the same value of k(X p ) phom — 1-7 
in eq.(107) gives perfect fits to the large a behavior 
of f L (a) for all 0.1 < r c < 0.8 (not shown). 

Figure ^ shows contours of the temperature of 
graphite dust in the ICM (T; cm ; top panel), in the 
clumps (T c ; middle panel) , and the temperature ratio, 
Ticm/Tc (bottom panel). The thick solid lines are the 



curves /l(ck) for which L c 



and the dashed 



lines are the curves fu (ex) for which M c = Mi cm - 
The contours of Ti cm wrap around /i,(a), whereas 
the contours of T c wrap around /m(ci). As a — > oo 
and / c — > 0, we find that T c decreases whereas Tj cm 
increases slightly. From the relation L oc MT 4+I3 , 
where f3 is the emissivity index, we get that 



T 



Tc 



M c 

Mir, 



L, 



(108) 
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Fig. 27. — Temperatures of graphite dust heated 
by uniformly distributed sources, as computed by the 
MGEP model over a wide range of clumpy medium 
parameters (ot,f c ). The top, middle, bottom panels 
show temperatures in the ICM, the clumps, and the 
ratio, respectively. The diamond and cross indicate 
results that were compared with MCRT. 



and with the conditions from eq.(104), this equation 
explains why the situation of /m(c>0 < fc < fh{ot) 
leads to T c < Ti cm . We can also explain why the 
temperature ratio must increase as a — > oo. Substi- 



tuting eq.(107) into eq.(103) gives 



M c 



(109) 



and where 



along the /l(o;) curve, where L c = 1 
C is a constant. Using this in eq. (|lbsj ) we get that 
along the /l (a) curve 



(110) 



The 



Therefore Ti cm /T c must increase as a — > oo. 
net result is that T c decreases faster than T icm in- 
creases because M = M c +M icm is a constant whereas 
L a bs = L c +L icm = 2L icm is decreasing as a — > oo and 
fc = /l(«). Note that the spectrum of the emitted 
radiation is dominated by dust at temperature T c if 
fc > /l(ck), since then L c > Li cm , or it is dominated 
by dust at temperature Ti cm if f c < /l(«). 

The temperature of silicates also exhibits the same 
variations but with lower values since silicates absorb 
only about 25% of what graphites absorb for the cho- 
sen composition. The same type of pattern, shown 
here for specific values of r c , T s , L s , and source type, 
is found for other values of r c , T s , L Sl and the cen- 
tral point source or uniformly illuminating external 
source types. The dust temperatures are lower as L s 
is decreased, and there is less variation in Ti cm /T c as 
T s or r c is decreased. Cases of T s > 20, 000K would 
require modeling the ionization of gas and the heating 
of dust by absorption of Lya photons. 

8. SUMMARY AND CONCLUSIONS 

We have introduced an analytical model for the es- 
cape and absorption of radiation in two-phase clumpy 
media, based on combining the mega-grains approx- 
imation of HP93 with escape/absorption probability 
formulae for homogeneous media, as summarized in 
§|], and referred to by the acronym MGEP. Enhance- 
ments of the mega-grains approximation developed 
in this paper include: (1) improved and extended the 
formula for effective optical depth to all clump fill- 
ing factors, (2) improved the formula for the effec- 
tive albedo, (3) developed a new approximation for 
the effective scattering asymmetry parameter, and (4) 
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developed new approximations for the fractions ab- 
sorbed by each phase of the medium. We also devel- 
oped a new approximation for the fraction of photons 
from a uniformly illuminating external source that are 
absorbed in a sphere of dust when scattering occurs. 

The space of parameters describing a two-phase 
clumpy dusty environment is six-dimensional: the 
clump filling factor, the clump to interclump medium 
(ICM) density ratio, and the clump radii, together 
characterize the morphology of the medium, and 
the interaction (absorption plus scattering) coeffi- 
cient, scattering albedo and asymmetry parameter 
arc the three optical parameters of the dust that 
vary with wavelength. The analytic MGEP model 
was compared with Monte Carlo simulations of ra- 
diative transfer (MCRT) over a subspace of the six- 
dimensional parameter space, and we found good 
agreement for most of the parameter values checked. 
More importantly, the qualitative behavior of the 
MGEP model agrees very well with MCRT simula- 
tions, giving us confidence in applying the MGEP 
model to parameter ranges that have not yet been 
checked by MCRT simulations. Three types of source 
distributions were studied: the central point source 
(C), uniformly distributed internal sources (U), and 
uniform illumination by external sources (X). Source 
type U is the least absorbed, whereas source type X 
is the most absorbed at low optical depths and source 
type C is the most absorbed at high optical depths. 

The MGEP model was shown to predict very well 
the emerging SED in a realistic simulation of starlikc 
sources heating a clumpy dusty medium, as compared 
to MCRT simulations. Furthermore, the MGEP 
method requires just seconds to compute a full spec- 
trum simulation, in comparison to hours of computa- 
tion when using the MCRT method. 

From MGEP simulations over a wide range param- 
eters characterizing the clumpy medium we find that 
for constant clump to ICM density ratio, a, the to- 
tal luminosity absorbed by the clumpy medium at- 
tains a minimum at the filling factor, f c , for which 
the luminosity absorbed by clumps (L c ) and the ICM 
(Licm) are equal. The curve of f c versus a for which 
L c — Li cm is found to be proportional to a - 2, a 
consequence of the clumps becoming optically thick, 
whereas the curve of f c for which the clumps and ICM 
have equal mass is proportional to a -1 , and these di- 
verging behaviors cause the temperature of dust in 
clumps to decrease as a — > 00 and f c — ► 0. Physically, 
the dust in opaque clumps shields itself from radia- 



tion, thus reaching a lower equilibrium temperature 
than dust in the ICM. The extra parameters gained by 
introducing dumpiness allows for modeling more un- 
usual relationships between the luminosity absorbed 
by dust and the resulting dust temperatures. 
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A. FINAL FLUX APPORTIONMENT IN 
MONTE CARLO 

Let N be the number of photons emitted as a 
group at the start of the Monte Carlo simulation and 
suppose after k iterations of traveling, escape, ab- 
sorption and scattering there are Nk photons remain- 
ing in the medium. Since after each scattering the 
flux weight of each photon is reduced by the albedo 
< co < 1, the actual flux remaining after fc scatter- 
ings is co k Nk- For clarity in the following derivation 
define the flux weight factor 

W k EE LU k . (Al) 

After one more iteration of traveling there are Nk+i 
photons remaining, each still having weight Wk, so 
that the flux escaping during the fc + 1 iteration is 
Ek+i = Wk (Nk — Nk+i). The remaining photons 
interact with the medium and are apportioned into 
absorbed and scattered fractions: the absorbed flux 
is Ak+i = (1 — of) Wk Nk+i, and the flux remaining 
after scattering is u> Wk Nk+i — Wk+iNk+i- Carrying 
this analysis forward by induction gives the general 
formula for the escaping flux after each iteration: 

Ek+i = Wk(N k ~Nk+i) 
Ek+2 = Wk+i{N k +i - N k+2 ) 

Ek+n = UJ^-^Nk+n-l-Nk+n); (A2) 

and the absorbed flux after each iteration: 
Ak+i = (1 - w) Wk Nk+i 

Ak+2 = (l-Cj)Wk+lN k +2 

Ak+n = (1 - w)w fc+n - 1 JV fc+n . (A3) 
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We seek expressions for the infinite sums of E k+n and 
A k+n over n > in order to terminate the Monte 
Carlo simulation at the fc-th iteration and still have 
an accurate estimate of the escaping flux. Assume 
that for all n > 



N, 



fe+Ti 



N, 







k+n-l 



(A4) 



so that after k iterations the fraction of photons re- 
maining after each additional iteration is assumed to 
have reached a steady state. Then we have 

N k+n = (3 n N k , (A5) 
JVfc+n-i - N k+n = (1 - PW^Nk . (A6) 

The sum of the escaping flux is then: 



n=l 



Ek+n 



fc+n-l 



N k+n ) 



oo 



n=l 

u k N k (l-l3)/(l~Lu[3) 



(A7) 



Since uj k N k is the known amount of flux remaining 
after k iterations, the fraction that would eventually 
escape if the iterations were continued is then 



f esc 

Jk — 



1-/3 



(A8) 



assuming that eq.(A4) is true. Similarly for the sum 
of the absorbed flux 



oo 

^ Ak+r, 
n=l 



(1 — Lj)uj 



N, 



cj k N k (l-Lj)J2t 



k-j-n 



u k N k (l-u)(3/{l-uf3), (A9) 



and it is clear that 



f esc i r abs 
Jk ' J k 



1-P , (l-w)/3 



I-lu/3 



1, (A10) 



giving a check of the derivations. Note that the sim- 
pler assumption of / fc esc = ui and f£ bs = 1 — ui corre- 
sponds to the case of /3 = 1/(1 + u>), which incorpo- 
rates no knowledge of the particular radiative transfer 
situation and could not possible work for all types of 



sources and geometries. In contrast, the final itera- 
tion escaping fraction given by eq. (|A8|) is based on a 
partial history of absorption and scattering and there- 
fore can better estimate the true value. 

To examine the effects of relaxing the assumption 
of eq.(A4), define the remaining fraction at the fc-th 
iteration: 

& = #- - ( A11 ) 

the ratio of the number of photons remaining after 
the k-th iteration over the number of photons in the 
medium before the k-th iteration. In our experience 
with Monte Carlo simulations we find that the se- 
quence (3 k (k — 1,2,3, ...) is monotonically either in- 
creasing or decreasing (does not oscillate), and the 
direction of convergence depends on the source geom- 
etry and the value of g, the phase function asymme- 
try parameter. In our Monte Carlo simulations with 
uniformly distributed internal emitters, the sequence 
Pk con verges rapidly, and using (3 k in place of (5 in 
eq.(A8) gives an f k esc that predicts the total escaping 
flux quite well. The remaining fraction j3 k is found 
to be an increasing sequence when g < g*(r), where 
g* (r) is approximately given by eq. (|2l]) , since photons 
tend to get trapped in the medium when scattering 
is more isotropic. When g > g*(r) then (3 k is a de- 
creasing sequence, since photons tend to escape after 
more scatterings if the angular scattering distribution 
is on average more forward. When g = g*{r) we find 
that the sequence f3 k is essentially constant. This de- 
pendence on the asymmetry parameter g is related 
to the validity of the es cape probability formula for 
scattering discussed in §3.2 and further analyzed in 
Appendix C2. The behavior of the sequence (3 k is 
essentially the same for a uniformly illuminating ex- 
ternal source of photons. In the extreme case of a 
single central source, it takes more scatterings for (3 k 
to converge and it is always a decreasing function of 
k. Thus / fe esc always slightly overestimates the escap- 
ing fraction for the case of a central source, but never 
underestimates it. When the scattering is more for- 
ward the rate of decreasing (3 k is faster during the first 
few scatterings, since in this case the likelihood that 
a photon will escape increases after each scattering 
more than when the scattering is isotropic. 

The discovery of the above equations and behav- 
ior was facilitated by the use of dynamic arrays to 
represent a group of photons, a standard feature of 
the Interactive Data Language (IDL), a product of 
Research Systems Incorporated (RSI). 
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B. 



RANDOM CLUMPS IN A TWO PHASE 
MEDIUM 



In section 2/2 we introduced equation (|llj) for the 
number of identical random clumps needed to pro- 
duce a given volume filling factor in a two-phase 
medium. Here we give a more detailed derivation of 
that equation. Let X be a randomly chosen point in 
the medium. Then the probability that a randomly 
placed clump will contain the point X is simply 



P ; 



Vc 

V 



(Bl) 



where v c is the volume of a clump, and V is the total 
volume of the medium. Continue to randomly place 
more clumps in the medium, without regard to over- 
laps, for a total of N c clumps. Then the probability, 
P(n), that the point X will be contained in n clumps 
is given by the binomial distribution: 



P(n) 



p n (l-p) 



N c -r 



(B2) 



In particular, the probability that X is not in any 
clump is 

P(0) - (l-p) N ° (B3) 

and so the fraction of the volume V occupied by the 
clumps is 



f c = 1-P(0) = 1 - (l-p) N * . 
Solving for N c gives the desired equation 

ln(l - f c ) 



In(l-f) 



(B4) 



(B5) 



for the number of identical randomly placed clumps 
that will have a total filling factor of f c , when each 
clump has a filling factor p. 

The expectation value for n of the binomial distri- 
bution [eq.dB^)] is 



E(n) = ^nP{n) = N c p 



V 



Qc , (B6) 



equaling the porosity of the clumps, Q c , which was 
introduced for the mega-grains approximation. As 
p — » and N c — ► oo the binomial distribution is well 
approximated by the P oisson distribution having the 



same expe ctation (e.g. |Bevington fc Robinson 1992 
bp. 23-28|) : 



P{n) 



(B7) 



For n = this approximation gives 
f c = 1 - P(0) ~ 1 - e 



(B8) 



providing another equation relating porosity and fill- 
ing factor, which is exact for our purposes, since 
p < 10~ 3 (e.g. when r c < 0.1 Rs ) and N c > 1. 
When f c -c 1 then of course f c ~ Q c . 

Another quantity of interest is the average number 
of clumps encountered along a randomly chosen line 
of sight, also called the covering factor, F c . Consider 
a cylinder of radius r c and length L centered on the 
line of sight. The volume of the intersection of the 
cylinder with the random clumps is irr^Lfc and we 
can estimate the number of clumps in the intersection 
by dividing by the volume of a clump: 



Fr 



4 s 
— ITT 

3 c 



3Jc 
4r c 



(B9) 



This is actually more like a lower bound estimate since 
the clumps may overlap. An upper bound is obtained 
by using the porosity in place of the filling factor: 



Fr. < L 



30c 



(BIO) 



C. ESCAPE AND INTERACTION PROB- 
ABILITIES FOR A SPHERE 

C.l. No Scattering 

Here we derive equations ( p6| ) and ([l9|), which are 
exact formulas for the interaction and escape proba- 
bilities of external and internal sources, respectively, 
in a homogeneous sphere of dust, ignoring the effects 
of scattering. Thus, scattering and absorption are 
considered together as interactions causing extinction 
of photons along a line of sight, and by extinction es- 
cape probability we mean the probability of escaping 
without being scattered or absorbed. Consider a ray 
intersecting the sphere at an angle 6 with respect to 
the surface normal vector, which we call the impact 
angle. The length of the chord created by the in- 
tersection is 2RcosO, where R is the radius of the 
sphere. Defining r = pnR, the optical radius, where 
K is the dust interaction coefficient and p is the den- 
sity of dust, the extinction optical depth of the chord 
is 2rcos0, and this will be used in the derivations. 

To derive the interaction probability, eq. fl26|) , for 
the case of a uniformly illuminating external source, 
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consider a beam of photons traveling in parallel rays 
impacting a hemisphere. Upon computing the trans- 
mission of the beam we can get the interacting frac- 
tion of photons, and then by symmetry this gives the 
interaction probability for all possible beams impact- 
ing the sphere, i.e. a uniformly illuminating external 
source. Let Iq be the intensity of each parallel ray 
in the impacting beam. Then the intensity emerg- 
ing from the other side is reduced by the extinction, 
which depends on the impact angle as follows: 



Iout(T,6) = I Q exp(-2rcos< 



(CI) 



The total flux emerging without interaction is com- 
puted by integrating over all impact angles with re- 
spect to solid angle: 



F out (T) 



2nl 



271-7, 



ttIq 



,tt/2 








/ e 


- 2rcose cos 


9 sin 6» 


dd 


Jo 
















/ A*e~ 


2T " dp 




(C2) 


Jo 








1 


(7+272-) 






272 _ 






5 



where we have used the substitution of variables 
p = cos9 and integration by parts. The flux that 
would emerge if the sphere was empty is computed 
by integrating with r = in eq.(|C2|), obtaining 



F = F out (0) 



2ttIq I pdp 
'0 



ttI 



So the fraction of photons that interact is 
F - F out {r) 



Pi(r) = 



1 



F 
1 

272 " 



(C3) 



(C4) 



1 

272 



giving the interaction probability for a uniformly illu- 
minating external source. For the case of an optically 
thin sphere note that the extinction behaves as 



1 — 2-r/i when r <C 1, 



(C5) 



and so the emerging flux is 
F o ut(r~0) : 



27rJo / p (1 — 2rp) dp 
Jo 



At 

¥ 



(C6) 



and as expected this gives 



4t 

P,(t~0) ~ — 



(C7) 



For the case of uniformly distributed internal emis- 
sion, the probability of escaping without interactions, 
eq.(pj|), can be derived in a similar fashion. Let e be 
the emission per unit volume per second. Then the 
non-interacting intensity emerging from a ray at an 
angle 9 with respect to the surface normal is 

Wr,0) = — (l-e~ 2TCOse ), (C8) 

obtained in the standard fashion by integrating the 
transfer equation with no scattering along the chord 
of optical length 2rcos6' through the sphere. The to- 
tal non-interacting flux emerging in any given direc- 
tion is computed by integrating with respect to solid 



angle, as in eq. 
F out {T) = 



2ire 

f)K 
27T£ 

pK 
7T£ 

pn 



over all the parallel rays: 

tt/2 



-It cos 6 



(1-e 

f (l-e- 2 ^) fidn 
Jo 



cos 8 sin 9 d9 



(C9) 



1 



1 

272 



2r 2 J 



-2t 



obtaining a result similar to the interaction proba- 
bility above, due to the same exponential term in the 
integral. If the medium had zero absorption and scat- 
tering, the intensity emerging from a ray is 



I o (R,0) = 2eRcost 



(CIO) 



The total flux emerging in any given direction is again 
computed by integrating with respect to solid angle 
over all the parallel rays: 



/>2"7T /* 7r /2 
F (R) = d<f> . 

Jo Jo 



= AneR / ifdpL = 



dd cos sin 9 1 (R, 9) 
AireR 



(Gil) 



So the fraction of photons that escape the sphere is 
F out {T) 



Pe(r) 



F (R) 
3 



(C12) 



4r 



1 



1 

272 



1 

272 



using the fact that r = pnR, giving the extinction 
escape probability for a uniformly distributed internal 
source. 
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Table C2 
Multiple Scattering and Escape 



# I interacting 



absorbed 



scattered 



escaping 



1 -P 

L0{1-P) 2 

u 2 (l - P) 3 



(1- W )(1-P) 
(1 - cu) cj (1 - P) 2 
{1~lu)lo 2 {1- P) 3 

(1- ^w"-^! -P) n 



w(l-P) 
w 2 (l-P) 2 
w 3 (l - P) 3 

w"(l- P)" 



Pw(l - P) 
Pw 2 (l - P) 2 

Pw 3 (l - P) 3 
Pu> n (l - P) n 



C.2. Including Scattering 



Equation (20), which approximately includes the 
effects of scattering into any extinction escape prob- 
ability for the case of uniformly distributed isotropic 
emitters in a bounded medium, was presented in Lucy 
et al. (1991) but no derivation was given. In this 
section we derive the equation and point out some 
interesting properties. The scattering albedo of the 
medium is the fraction of the extinction optical depth 
that is due to scattering: ui — T S cat/T e xt, where 
Text — Tabs + Tscat, is the sum of absorption and scat- 
tering optical depths from the center to the boundary 
of the medium, e.g. the radius of a sphere. Assume 
that we are given the escape probability, P e (r) , for ab- 
sorption only as a function of the optical depth of the 
bounded medium. We can immediately use P e (j) to 
give the fraction of emitted photons that escape with- 
out interacting (not absorbed or scattered) , by using 
t = r ext . For convenience we shall identify P = P e (r) 
in the following derivation. Thus a fraction P of pho- 
tons escape with no interactions and a fraction 1 — P 
interacts with the medium. Then by definition of the 
albedo, a fraction (1 — w)(l — P) is absorbed and the 
fraction w(l — P) is scattered for the first time. As- 
sume that the scattered photons are distributed uni- 
formly in the medium and the directions of travel are 
isotropic. Then we can apply the extinction escape 
probability again to find that a fraction Pw(l — P) 
escapes after one scattering. Repeating these steps, 
observe that a fraction w(l — P) 2 interacts with the 
medium to result in a fraction u> 2 (l — P) 2 that scatters 
a second time. Continuing to assume that the scat- 
tered photons are uniformly distributed and isotropic 
it follows that a fraction Pw 2 (l — P) 2 escapes after 
two scatterings. 



Table |C2| summarizes the analysis of the fractions 
in each state and extends it by induction. The first 
column is the number of interactions that have oc- 
curred. The column labeled "scattered" feeds back 
into the "interacting" column of the next row with an- 
other factor of 1 — P to give the next interaction. Note 
that for the first interaction, the photons that are 
absorbed have experienced zero scatterings, so that 
the number of scatterings for fractions in the "inter- 
acting" and "absorbed" columns is one less than the 
interaction number. The induction presented in the 
table shows that a fraction Pw"(l — P) n escapes af- 
ter n scatterings. Summing these escaping fractions 
over an infinite number of scatterings gives the sought 
after formula for the total escaping fraction: 

oo 

V esc (r,0j) = P^>"(l-P)" 



n=0 



P 



1 - w(l- P) 



(C13) 



Recall that the formula is valid (most accurate) at a 
single value of the scattering asymmetry parameter, 
which in the case of spherical geometry is approxi- 
mately given by eq. ( |2l| ) . The absorbed fractions listed 
in the third column of Table C2 can also be summed 



to get the total absorbed fraction: 



V abs {r,oj) = (l-w)(l-P)j^w n (l-P) T 



(l-co)(l-P) 
l-w(l-P) 



(C14) 



finding that P esc (r, cj) +V a b s (T,u>) = 1, which gives 
a check of the above derivations. 
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Let us define distributions, p e (n) and p a (n), for the 
probability that a photon will escape or get absorbed 
after n scatterings. To arrive at a probability distri- 
bution, the fraction Pui n (l — P) n that escapes after 
n scatterings must be normalized by dividing by the 
total escaping fraction, obtaining 



Pe(n) 



Puj n {l-P) n 

V esc (T,Uj) 

= [l-w(l- P)}u n (l- P) n , (C15) 



after substituting eq.( p!3 ). As mentioned before, the 
number of scatterings occurring for the absorbed frac- 



tions listed in the third column of Table C2 is actually 
n — 1, where n is the number of interactions given in 
the first column. Therefore, the probability distribu- 
tion of absorption after n scatterings is 



Pa(n) 



(1 -cj)w n (l -P) 

Tabs (T,U>) 

[l-o;(l-P)]o; n (l 



n+l 



P) n , (C16) 



resulting in the fact that p a (n) — p e (n) for all n. 
Of course this can be true only when P esc (r, cj) is 
valid. Monte Carlo simulations verify that the escap- 
ing and absorbed probability distributions are equal 



only when equation (C13) agrees with the Monte 
Carlo results, which occurs only for a single value of 
the scattering asymmetry parameter g, as discussed 
in section 3.2 for the case of spherical geometry. 

Defining A = 1 — P, the average number of scat- 
terings that a photon experiences before escape or 
absorption is calculated as 



(n sca t) 



^2,np e {n) 

n=0 

oo / oo k 

E £*(»)-£ 



Pe(n) 



k=0 \n=0 



1 - LU k+1 A k+1 



1-uA 



X> fc (l-P)* 



fc=i 



1 



1 - w(l- P) 
V esc {r,u) 



Pe{r) 



- 1 



(C17) 



Monte Carlo simulations for the case of uniform emis- 
sion within a homogeneous sphere of absorbers and 



scatterers again verifies that when the asymmetry 
parameter takes on the particular value g — <?*(t), 
approximately given by eq.(^l|), then the average 
number of scatterings experienced by photons is the 
same whether the final state is escape or absorption: 
(n S cat) esc = (nscat}* = (n SC at) abs - For other values 
of g we find that (n scat ) esc < (n sca t)* < (n SC at} abs 
when g < g* (r) , so that absorption occurs after more 
scatterings than escape since there is more absorp- 
tion than predicted by eq.(C13), whereas (n sca t) esc > 
(nscat)* > (n S cat) abs when g > g*(r), since there are 
more escaping photons than predicted. We conjec- 
ture that the ratio of ( 

nscat) eS c ^ (nscat) abs he 
affected only by g and r because the ratio is depen- 
dent only on geometry: the directions of scattering 
and proximity to the boundary of the medium. The 
albedo u> gives probability of scattering relative to 
absorption and so affects the magnitude of (n SC at)*, 
but it does not enter into the geometry of photon 
paths, possibly explaining why g*(r) is independent 
of u>. Note that as the medium becomes optically 
thick (t — * oo) then P — * and then eq.(C17) pre- 
dicts that (nscat)* — * (1 — ui) from below, and so for 
.9 < 9*( T ) w e have that (nscat) esc < w/(l — u>), which 
is effectively for all g since g*(r) — > 1 as r — > oo. 



D. ANALYSIS OF CLUMP OVERLAPS 

To extend the mega-grains approximation to high 
filling factors, consider that as f c — + 1 it is the increase 
of clump overlaps that makes the mega-grains model 
become unrealistic. Since we ignore the extra density 
occurring in overlaps, the effective radii of the clumps 
can be reduced to eliminate most of the overlapping 
and the number of clumps then must be increased to 
retain the same filling factor. The volume of clump 
overlaps is V (Q c — f c ) and so the average overlap per 
clump is 



v(Qc-.fc) 



-fc 



(Dl) 



where we have applied eq. (|44|) and v c is the volume of 
just one clump. Suppose there is a pair of clumps that 



overlap in volume by the amount given in eq.(Dl), 
then to find the reduced radius at which the pair of 
clumps would not overlap we must solve 



I c {' Qc - fc 

2 



Qc 



h) (D2) 
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for h, the half-depth of the overlap, where r c is the 
clump radius. The reduced radius that eliminates the 
overlap is then given by r c —h. We have solved eq. ( |E>2] ) 
numerically for the full range of filling factors < 
f c <l and find that r c (1— f c ) < r c — h as f c — > 1, thus 
eliminating the overlap. We propose to substitute 
r c (l — fc) 1 for all instances of r c in the mega-grains 
equations, where 7 is an optional tuning parameter. 
In this new model, when 7 = 1, the density of clumps 
is 

Hc ~ 4^(1 - / c ) 3 Airrl ( ' 

thus greater than originally defined, and we assume 
that the radius renormalization has effectively elimi- 
nated overlaps while preserving the filling factor. 



E. DISTRIBUTION OF DUST TEMPERA- 
TURES AROUND A POINT SOURCE 
OF RADIATION 

Given an isotropic point source of radiation in a 
spherically symmetric medium, the resulting distri- 
bution of radiation and heating of the dust is then a 
function of only the radial distance r from the point 
source. We will show that the probability distribution 
of dust temperatures can be approximated by a power 
law function if the dust density and absorbed lumi- 
nosity versus r, and if the emitted luminosity versus 
dust temperature can be approximated by power law 
functions. Letting N{r) be the number of dust grains 
in a sphere of radius r, our approach is to obtain ap- 
proximations for 



dN _ dN (dT 
dT dr \ dr 



(El) 



which then gives the functional form for the distribu- 
tion of temperatures. 

Assume that the luminosity absorbed by the dust 
at distance r varies like an inverse power law with 
exponent 77, 



L a bs{r) ~ m(r) r n , 



(E2) 



which in the optically thin case is nearly exact with 
7] = 2, where m{r) is the mass of dust in a thin shell at 
radius r. When the dust is optically thick we expect 
r) > 2, because the absorbed luminosity decays more 
rapidly than any power law (77 is then also a function 
of r) . The luminosity emitted by the dust scales with 



the dust temperature, T(r), approximately as 

/>oo 

L em [T(r)} = 47rm(r) / B v [T{r))K v dv 
Jo 

~ m(r) T(r) 4+/3 , 



(E3) 



if the emissivity per unit mass of the dust scales with 
frequency like k„ ~ zV 3 , and usually < (3 < 2. 
Equating the absorbed and emitted luminosities of 
the dust we get an approximate relationship between 
dust temperature and radial distance: 



T(r) ~ r 



(E4) 



The rate of change in temperature with respect to 
radial distance then varies as: 



dT 
dr 



1+7 



T 



(E5) 



where we have also inverted eq.(E4) and substituted 
for r. 

Assume that the dust density is approximately a 
power law p(r) ~ r~ 5 , so that the number of dust 
grains, N(r), in a sphere of radius r also varies like a 
power law function: 



N(r) = 4tt / p(s) s 2 ds 



„3-s 



Then 



dN 
dr 



„2-<5 



T 



(E6) 
(E7) 



where we have again substituted for r using eq.(E4). 
Combining eqs.(E5) and (E7) into (El): 



^1 „ T^-M^ri- 1 -^) 



dT 



T C*-3)(±±fi)-l 



(E8) 



From this scaling approximation we presume that the 
distribution of dust temperatures follows a power law 



p(T) = aT» 



with exponent 



=(,i-:!,('i±^|-l 



(E9) 



(E10) 



which is certainly negative as long as S < 3. The 
normalizing factor a is determined by requiring 



p(T)dT 



1 



(Ell) 
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where T max is usually the dust sublimation tempera- 
ture (~ 1500K) and T m j n is the minimum dust tem- 
perature. Actually, the minimum dust temperature 
is the major free parameter in this theory, and is de- 
termined by balancing the absorbed and emitted dust 
luminosities: 

r R /-T max 

L abs {r)r 2 dr = M p{T)K{T)dT, (E12) 

J ^min 

where M is the total mass of dust contained in a 
sphere of radius R, and where 

K(T)= / B v {T)n v dv (E13) 
Jo 

is the Planck averaged emissivity at temperature T. 
Most of the emission comes from dust at temperatures 
near T m i„ since those dust grains occupy most of the 
volume and the dust temperature decays rapidly with 
distance from the point source. Once r m ; n is deter- 
mined, the IR emission spectrum from the dust grains 
is given by 



F v = 4itMk u 



B v (T)p(T) dT 



(E14) 



The behavior of the emissivity of graphite or sili- 
cates changes at wavelengths longer than about 10/im, 
therefore we find that the power law approximation 
of the emissivity averaged with a Planck function at 
a given temperature (proportional to the luminosity 
emitted by dust), 



K(T) ~ T 4+ft , 



(E15) 



requires two exponents, Pi, with i = 1 for T > Tb, 
and i = 2 for T < Tb, where Xj is called the emissivity 
break temperature, corresponding to the wavelength 
at which the behavior of the dust emissivity changes. 
From examination of the behavior of K(T), the values 
of Pi and T b are found to be approximately 



and 



(1.0, 2.0) (graphite) 
(0.5, 2.0) (silicates) 



80 K (graphite) 
150K (silicates) 



(E16) 



(E17) 



Consequently, the power law distribution of temper- 
atures around a central source will have an exponent 
depending on the dust temperature range: 



Pi 



(5-3) 



4 + Pi 



- 1 



(E18) 



yielding power law probability distributions 

Pi (T) = a{T^ , (E19) 

where the constants aj are determined by requiring 
thatpi(T b ) = P2 (T b ) and 



/ P i{T)dT + p 2 (T)dT=l, (E20) 

J Tb J T m i n 



yielding (defining ^ = p { + 1) 



a 2 = ai T b 7: 



(E21) 



ai 



7i 72 



72^ + (7i - 72 W 1 + 7i^ n r 6 71 " 72 



This approximation for the temperature distribution 
and resulting IR emission spectrum were found to be 
in reasonable agreement with Monte Carlo simula- 
tions for the case of S = 0. 
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